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13  ABSTRACT  ( Maximum  200  words) 

A  seedling  program  was  conducted  to  determine  whether  load  information  is  stored  in  the  fracture  surface 
features  of  failed  aircraft  components  and,  if  so,  to  explore  methods  for  extracting  this  information. 

Fracture  surfaces  produced  under  systematically  varied  cyclic  load  conditions  in  laboratory  specimens  of 
titanium  turbine  blade  alloy  were  provided  to  the  program  by  an  aircraft  engine  manufacturer.  Three- 
dimensional  topographs  of  the  surfaces  were  obtained  via  confocal  optics  and  miscroscopy  and  analyzed 
with  Founer  and  wavelet  techniques.  Both  analysis  techniques  confirmed  the  existence  of  load 
information  in  the  fracture  surface  topography.  Moreover,  correlation  was  found  between  Fourier  results 
and  fatigue  load  parameters. 

The  topographic  characterization  and  analysis  method  was  then  applied  to  an  F-118  engine  component 
that  had  cracked  in  service  to  determine  whether  the  load  history  that  caused  the  cracking  could  be 
extracted.  Using  the  results  from  the  laboratory  specimens  as  a  reference  database,  the  method  provided 
a  quantitative  estimate  of  an  influential  fatigue  load  parameter  and  indicated  its  variation  with  crack  growth 
distance. 
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13.  ABSTRACT  ( continued ) 

The  results  from  a  companion  component  were  in  agreement.  Thus,  a  route  has  been  demonstrated  for 
extracting  load  spectra  information  from  fracture  surfaces,  and  away  is  pointed  for  a  leapfrog  advance  in 
the  state  of  practice  of  failure  analysis. 

A  plan  is  presented  for  elevating  fractography  from  a  qualitative,  subjective,  two-dimensional  art  to  a 
quantitative,  objective,  three-dimensional  science.  The  resulting  methodology  will  enable  load  conditions 
and  crack  evolution  kinetics  to  be  extracted  from  the  fracture  surfaces  of  a  failed  part,  thereby  providing 
aircraft  component  designers  with  previously  unobtainable  critical  design  information  that  will  enhance  air¬ 
readiness  and  avoid  hundreds  of  millions  in  annual  maintenance  costs. 


iii 


SUMMARY 


Cracks  that  form  in  aircraft  engine  turbine  blades,  disks,  vanes,  and  rotors  under  high 
frequency  vibrational  loads  during  flights  have  caused  many  Class  A  incidents  (loss  of  aircraft)  in 
the  USAF  fighter  fleet.  These  high-cycle  fatigue  failures  are  responsible  for  unacceptable  loss  of 
life,  costly  redesigns  and  retrofits,  and  huge  maintenance  costs.  More  important,  high-cycle 
fatigue  failures  result  in  grounding  of  the  fleet  for  days  to  weeks.  The  USAF  has  identified  high- 
cycle  fatigue  as  the  primary  air  readiness  problem. 

The  high-cycle  fatigue  problem  exists  because  the  load  spectrum  experienced  by  engine 
components  is  not  known.  In  the  damage  tolerance  philosophy  used  to  manage  aircraft  engine 
integrity,  this  means  crack  growth  rates  cannot  be  adequately  predicted,  which  means  components 
cannot  be  reliably  designed  and  inspection  intervals  cannot  be  set  with  confidence. 

Thus,  a  way  is  needed  to  determine  the  load  history  of  in-service  components.  The  load 
spectra  deduced  by  finite  element  analysis  of  engine  vibration  characteristics  are  not  pertinent 
because  of  aberrations  from  mission  uniqueness,  unforeseen  pilot  maneuvers,  wind  gusts,  and  the 
like.  Pilot  reports  are  subjective  and  on-board  sensors  are  not  usually  fielded.  A  possible  solution 
is  to  extract  load  spectra  information  from  the  markings  on  the  crack  surfaces.  This  seedling 
program  investigated  whether  parameters  of  the  fatigue  load  could  be  obtained  by  analyzing  the 
crack  face  topography. 

We  obtained  laboratory  specimens  of  Ti-6AMV  that  had  been  tested  to  failure  under 
prescribed  fatigue  loading  conditions  and  selected  areas  of  the  fracture  surfaces  that  had  been 
produced  at  various  values  of  Knm,  stress  range,  and  stress  ratio — the  three  important  load 
parameters  for  fatigue.  Using  a  confocal  optics  microscope,  we  made  topographic  maps  of  these 
areas  and  analyzed  the  topographs  with  Fourier  and  wavelet  techniques. 

The  energy  power  spectral  density  (EPSD)  curves  from  the  Fourier  analysis  showed  clear 
systematic  variations  with  AK,  stress  ratio,  and  hence  K,^.  The  curves  varied  in  height  with  AK; 
the  shape  of  the  curves  varied  with  stress  ratio.  For  any  given  wavelength  between  0.5  and  10  pm, 
a  linear  relationship  existed  between  the  square  root  of  EPSD  and  AK.  Moreover,  a  stress  ratio  of 
0.8  produced  a  linear  EPSD  curve  in  this  wavelength  range,  whereas  a  stress  ratio  of  0.1  produced 
a  concave  shape. 
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These  results  show  that  load  parameter  information  is  encrypted  in  the  roughness  details  of 
the  fracture  surfaces.  Wavelet  results  confirmed  these  findings,  but  were  less  definitive  because  the 
fracture  areas  examined  were  small  relative  to  the  grain  size  of  the  alloy  and  because  of  noise  in  the 
elevation  data.  Nevertheless,  the  Fourier  results  reveal  the  opportunity,  when  fully  explored  and 
properly  analyzed,  to  determine  quantitative  load  spectra  parameters  from  the  fracture  surface 
topography. 

We  then  examined  an  aircraft  component  that  had  acquired  a  crack  during  service  to 
determine  if  the  load  history  that  caused  the  crack  could  be  identified.  The  topographies  of  areas 
along  the  crack  growth  direction  were  characterized,  then  analyzed  with  the  Fourier  method.  The 
EPSD  curves  from  the  components  were  interpreted,  using  the  results  from  the  laboratory 
specimens  to  estimate  quantitatively  the  load  conditions  responsible  for  the  failure  and  their 
variation  with  crack  depth.  Strictly  from  the  fracture  surface  elevation  data,  we  determined  that  the 
crack  grew  to  a  length  of  1 .5  mm  under  rather  constant  AK  level  of  about  5  MPaVm  after 
initiating  at  a  higher  value.  At  this  crack  depth,  the  AK  increased  strongly  and  linearly  to 
10  MPaVm  in  the  next  millimeter  of  crack  advance.  We  further  estimated  that  the  stress  ratio  was 
about  0.8. 

Thus,  a  route  has  been  demonstrated  for  extracting  load  information  from  fracture  surfaces. 
The  next  steps  are  to  develop  topography  characterization  techniques  capable  of  resolving  finer 
fracture  surface  features,  develop  size  scaling  procedures  for  roughness  features  ranging  from 
nanometers  to  millimeters,  develop  wavelet  analysis  techniques  to  obtain  fracture-site-specific 
signals,  and  develop  a  computational  model  of  microstructural  deformation  and  microfailure.  The 
resulting  technology  will  enable  load  conditions  and  crack  evolution  kinetics  to  be  extracted  from 
the  fracture  surfaces  of  a  failed  part  and  will  constitute  a  leapfrog  advance  in  the  state  of  practice  of 
failure  analysis. 
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INTRODUCTION 


Fractography.  the  art  of  interpreting  markings  on  fracture  surfaces  to  deduce  how  and  why 
a  failure  occurred,  is  an  important  aspect  of  many  failure  investigations.1  Fractography  is  often 
relied  on  to  determine  whether  the  material  was  defective,  the  component  design  was  inadequate, 
or  the  allowable  operating  conditions  were  exceeded.  Hence,  fractography  implicates  wha*.  was 
responsible  for  the  failure  and  provides  the  basis  for  modifications  to  reduce  the  likelihood  of 
future  failures. 

Generally,  fractographic  results  are  qualitative  and  consist  of  the  failure  mode,  the  load 
conditions,  the  fracture  initiation  site,  and  the  propagation  direction.  Fractures  are  classified  as 
ductile,  brittle,  or  intergranular;  the  loads  that  caused  the  failure  are  identified  as  cyclic,  overload, 
or  environment-assisted;  the  crack  is  determined  to  have  initiated  at  a  material  defect,  a  machining 
or  damage  mark,  or  at  a  stress  concentrating  feature  of  the  component  geometry.  When  the  origin 
of  the  crack  can  be  pinpointed,  the  initiation  site  can  be  examined  with  microscopes  and  chemical 
analysis  instruments  to  determine  the  nature  of  the  defect. 

Although  all  this  information  is  very  useful,  the  failure  analyst  would  often  like  to  know 
additional  and.  if  possible,  quantitative  information.  For  example,  what  were  the  loads  (their 
magnitude,  duration,  frequency)  that  operated  and  led  to  the  fracture?  How  long  did  it  take  for  the 
crack  to  form?  And  once  formed,  what  was  its  growth  history?  Did  it  propagate  at  a  uniform 
speed,  or  did  it  experience  periods  of  acceleration,  deceleration,  or  perhaps  hesitation?  A  way  to 
quantitatively  extract  the  load  parameters  and  the  crack  kinetics  from  the  fracture  surfaces  would  be 
of  great  utility. 

Several  years  ago.  we  began  thinking  whether  such  information  might  be  encrypted  in  the 
fracture  surfaces  and.  if  so,  how  it  might  be  extracted.  We  participated  in  a  Navy  workshop2 
where  the  topic  of  interest  was  fatigue  failures  of  helicopters  and  fixed-wing  aircraft. 
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An  unacceptable  frequency  of  accidents  was  causing  safety  concerns,  incurring  huge 
maintenance  costs,  and  compromising  the  air  readiness  of  the  fleets.  Component  designers  at  the 
workshop  told  us  that  fatigue  problems  could  be  reduced  if  the  fatigue  load  spectra  were  known. 
Components  were  designed  on  load  spectra  deduced  by  finite  element  analysis  of  engine  vibration 
characteristics.  In  reality,  however,  the  load  spectra  experienced  by  a  component  during  service 
can  be  quite  different  because  of  mission  uniqueness,  unforeseen  pilot  maneuvers,  wind  gusts,  and 
the  like.  Overloads  are  not  reliably  obtained  from  pilot  reports,  and  on-board  sensors  are  not 
usually  fielded. 

It  occurred  to  us  at  this  workshop  that  load  spectra  might  be  obtainable  from  the  fracture 
surfaces.  We  reasoned  that  the  important  loads  in  the  spectrum  were  those  that  affected  the  fatigue 
process,  and  since  these  loads  affected  crack  growth,  they  might  manifest  themselves  on  the 
fracture  surfaces.  We  considered  that  the  extent  and  nature  of  out-of-plane  deformation  should  be 
a  function  of  the  loads  and  hence  should  be  accessible  from  the  topography.  Thus  the  fracture 
surface  may  act  as  a  load  spectrum  gauge! 


It  is  difficult  to  estimate  the  total  annual  cost  of  aircraft  fatigue  failure  to  the  military.  When  we  asked  the  question 
to  representatives  of  the  military  and  industry,  we  learned  that  the  costs  fall  into  several  categories,  failure 
investigations,  redesign  and  retrofitting,  and  maintenance,  but  the  greatest  concern  of  aircraft  fatigue  is  air  readiness. 
A  dollar  value  cannot  be  placed  on  having  the  fleet  grounded. 

Failure  investigations.  “The  costs  of  engineering  root  cause  investigations  frequently  run  into  the  multiple 
millions  of  dollars.  These  investigations  require  large  cross-disciplinary  groups  of  investigators  (from  20  to  100 
people),  laboratory  verification  of  failure  scenarios,  and  costly  engine  tests  (about  $40,000  per  day  ).  (One  specific 
accident  investigation  cost  $23M.  Otha  Davenport,  SPO  USAF.)  These  investigations  also  can  take  many  months 
to  complete,  during  which  time  the  unresolved  problem  can  lead  to  additional  field  failures,  possibly  resulting  in 
loss  of  aircraft  and  life."  Dr.  Jerrol  W.  Littles,  Jr.,  Pratt  &  Whitney 

Redesign  and  retrofitting.  “The  potential  for  saving  a  costly  redesign  and  retrofitting  of  the  fleet  if  the  cause  were 
known  and  another  solution  such  as  damping  could  be  used  would  be  enormous.  For  example,  the  fourth  stage 
blade  failure  in  a  single  crystal  material,  which  was  never  fully  understood,  ended  up  costing  the  government 
probably  over  20  million  dollars.  If  we  knew  the  failure  mode,  we  might  have  prevented  failure  by  adding  damping 
or  redirecting  air  flow ."  Dr.  Theodore  Nicholas,  AFRL,  Program  Manager,  Materials  Characterization,  HCF 
Science  and  Technology  Program.  The  Navy  experienced  “extreme  trouble  w  ith  cracking  in  helicopter  rotor  heads 
caused  by  uncertainty  in  fatigue  crack  growth  rates,  because  the  load  spectrum  was  not  known.  The  problem  was 
addressed  at  great  expense  by  redesigning  the  component  and  changing  the  material.  John  T.  Cammett,  Naval 
Aviation  Depot,  Cherry  Point,  NC. 

Maintenance  costs.  Annual  costs  for  just  one  specific  failure  mode  and  component — high-cycle  fatigue  of  aircraft 
engines— were  $200M  in  FY94,  $150M  in  FY92-3,  and  S220M  in  FY88-91.  Robert  J.  May,  Jr.,  USAF 
Propulsion  Product  Group  Manager.  The  Navy,  with  a  fleet  of  10(Vs  of  helicopters,  requires  each  to  be  inspected  at 
intervals  as  short  as  every  10  hours  of  flight.  John  T.  Cammett,  Naval  Aviation  Depot,  Cherry  Point,  NC. 


Air  readiness.  High-cycle  fatigue  is  the  No.  1  air  readiness  concern  of  the  USAF.  Dr.  Deutch  USAF.  If  just  one 
recurring  failure  mode  could  be  eliminated,  countless  days  or  weeks  of  not  being  in  a  readiness  posture  in  the  field 
could  be  avoided. 
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This  idea  may  not  be  so  far-fetched.  For  example,  in  certain  situations  a  procedure  for 
deducing  AK  from  fracture  surfaces  already  exists.34  When  striations  are  visible  and  a  correlation 
with  load  cycle  and  striation  can  be  established,  the  rate  of  change  of  crack  length  with  load 
application  can  be  obtained  and  correlated  with  data  from  controlled  fatigue  crack  growth  tests  in 
the  laboratory  to  deduce  AK.  Even  in  the  absence  of  striations,  skilled  fractographers  can  often 
estimate  the  crack  growth  rate  to  an  accuracy  of  half  an  order  of  magnitude-  and,  using  the  Pans 
law,  obtain  an  estimate  of  AK.  Furthermore,  recent  work  by  fracture  scientists  in  Japan6 ' 
suggests  a  relationship  between  striation  height-to-width  ratio  and  the  stress  ratio,  R. 

To  be  more  useful  and  generally  applicable,  however,  the  load  spectra  extraction  procedure 
needs  to  be  more  systematic,  more  quantitative,  and  more  precise.  The  fact  that  engineers 
experienced  at  assessing  fracture  surfaces  can  reach  consensus  on  approximate  load  conditions 
suggests  that  an  objective  method  is  possible.  If  the  eye  and  brain  can  perform  this  function,  it 
must  be  possible  to  develop  an  instrument  and  software  to  do  the  same.  A  preliminary  effort  was 
conducted  by  SRI  to  test  whether  fracture  surface  topography  was  sensitive  to  load  conditions,8-9 
and  a  seedling  program  was  initiated  by  DARPA  to  conceive  an  approach  and  investigate  this 
possibility.  This  report  describes  the  results  of  the  seedling  program  and  presents  a  plan  for 
elevating  fractography  from  a  qualitative,  subjective,  two-dimensional  art  to  a  quantitative, 
objective,  three-dimensional  science. 
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BACKGROUND 


We  first  consider  how  a  fracture  surface  forms.  Figure  1(a)  shows  the  slip  bands  that 
formed  in  the  grains  of  a  nickel  alloy  after  1200  stress  cycles.10  These  bands  of  localized 
dislocations  lie  on  crystallographic  planes  of  low  shear  resistance  and  hence  are  oriented  differently 
in  adjacent  grains.  Where  the  bands  encounter  a  grain  boundary',  dislocation  motion  stops, 
because  of  the  misalignment  of  the  slip  planes  in  the  adjacent  grain.  Stresses  at  such  intersections 
are  created  because  the  displacement  produced  by  the  band  is  not  transferred  to  the  neighboring 
grain  and  the  strain  discontinuity  is  accommodated  in  the  grain  boundary. 

Figure  1(b)  shows  the  same  area  after  4000  cycles.  The  density  of  the  slip  bands  has 
increased,  and  the  resulting  increase  in  stress  at  the  grain  boundaries  has  caused  the  boundaries  to 
separate.  The  failed  grain  boundaries  become  islands  on  the  fracture  surface.  Moreover,  the  grain 
boundary  facets  may  contain  steps  where  the  slip  bands  intersected. 

As  cycling  continues,  slip  band  activity  continues,  additional  grain  boundaries  fail,  and 
ligaments  between  grains  stretch,  neck,  and  rupture,  resulting  in  a  macrocrack  and  two  opposing 
fracture  surfaces  containing  roughness  features  on  a  range  of  size  scales.  The  tortuous  topography 
belies  the  deformation  and  failure  processes  that  took  place  in  the  microstructure.  Second  phase 
particles,  inclusions,  pores,  and  inherent  flaws  (when  they  are  present)  also  influence  the  failure 
process  and  hence  the  fracture  surface. 

In  the  case  of  fatigue  failure,  an  additional  phenomenon  can  affect  the  fracture  surface.  In 
the  unloading  phase  of  the  load  cycle,  the  crack  faces  relax  toward  each  other.  If  they  make 
contact,  they  can  modify  the  roughness  features.  Especially  for  conditions  where  the  lower  loads 
are  a  small  fraction  of  the  maximum  loads  or  even  compressive,  the  features  standing  up  in  relief 
can  be  deformed  to  different  shapes.  This  phenomenon  is  called  closure  and  has  an  important 
effect  on  fatigue  crack  growth  rate. 

Thus,  the  topography  of  a  fracture  surface  is  dictated  by  the  microstructure  and  the 
deformation  and  microfailure  processes  exacted  on  the  microstructure  by  the  applied  loads.  This 
suggests  that  the  load  history  experienced  by  a  failed  component  is  recorded  in  the  fracture  surface 
and  encrypted  in  the  roughness  details.  The  challenge  is  to  develop  a  procedure  for  extracting  this 
information. 
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Figure  1 .  Fatigue  cracks  in  nickel  after  1200  (< 
load  cycles.  (Ref.  1) 
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NM-31 2581-1 


and  4000  (b)  tension-compression 


The  first  steps  were  taken  in  this  seedling  program.  To  determine  whether  fracture  surfaces 
were  sensitive  to  loading  conditions,  we  obtained  fracture  specimens  tested  under  cyclic  loads 
under  controlled  conditions  of  different  AK,  quantified  the  topographies  of  their  fracture  surfaces, 
and  applied  Fourier  and  wavelet  analysis  techniques  to  the  elevation  data.  We  looked  for 
correlations  in  roughness  signatures  and  features  of  the  analyses  with  load  parameters.  To 
demonstrate  a  route  for  extracting  load  parameters  from  fracture  surfaces,  we  then  performed  a 
similar  analysis  of  a  field  failure  and  applied  the  modest  database  from  the  laboratory  specimens  to 
estimate  AK  as  a  function  of  crack  length. 
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ANALYSIS  OF  LABORATORY  FRACTURE  SPECIMENS 


Pratt  &  Whitney,  which  has  conducted  many  fracture  tests  over  the  years  to  determine  the 
fatigue  properties  of  turbine  engine  material,  made  its  large  inventory  of  tested  specimens  available 
to  this  program.  We  selected  for  analysis  two  compact-tension  specimens  of  Ti-6A1-4V  that  had 
been  tested  to  examine  fatigue  crack  growth  behavior  in  the  near-threshold  regime.  A  crack  had 
been  driven  through  the  specimens  under  a  range  of  well-controlled  and  documented  loading 
conditions  according  to  the  ASTM  E647  requirements  for  load  shedding.11  The  fracture  surfaces 
were  produced  at  constant  stress  ratio,  R  =  K^/K^,  and  varying  stress  range,  AK  = 

-  KjHjn-  with  successive  areas  produced  under  a  slightly  lower  or  slightly  higher  AK.  Two 
values  of  R  were  examined:  0.1  and  0.8.  The  variation  of  load  conditions  for  each  specimen  is 
shown  in  Figure  2.  In  total,  13  load  conditions  and  1 19  fracture  surface  areas  were  examined. 

SURFACE  TOPOGRAPHY  CHARACTERIZATION 

Small  areas  on  the  fracture  surface  in  each  fatigue  test  section  were  selected  at  a  certain 
spacing  along  the  crack  growth  direction,  and  their  topography  was  characterized  by  a  system 
known  as  a  FRASTAscope.121'  This  system  was  developed  at  SRI  and  consists  of  a  confocal- 
optics-based  scanning  laser  microscope:  a  precision,  computer  controlled  x-y-0  stage;  and 
computer  programs  that  control  data  acquisition,  manipulation,  and  display.  Tables  1  and  2 
summarize  the  locations  of  the  areas  we  chose  for  examination  and  the  corresponding  AK  values. 
The  size  of  each  area  was  140.3  pm  in  the  crack  growth  direction  and  466.5  pm  in  the  orthogonal 
direction.  The  number  of  elevation  data  points  in  the  crack  propagation  direction  was  600  and  in 
the  orthogonal  direction  2000.  Data  spacing,  thus,  was  0.233  pm  in  both  directions.  Figure  3 
shows  the  contrast  image  and  gray-scale  elevation  image  of  a  typical  area. 
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Figure  2.  Stress  intensity  range,  AK,  as  a  function  of  crack  length  during  multiple 
fatigue  tests  in  Specimens  8353  and  8354. 


Table  1.  Crack  Length  and  Stress  Intensity  Range  for  Several 
Fatigue  Tests  in  Specimen  No.  8353  (R  =  0.1) 


File 

No. 

Test  Run 
Number 

Crack  Length 
(mm) 

Stress  Intensity 
Ranqe,  AK,  (MPaVm) 

1 

Run  No.  1 

9.9238 

8.52 

2 

10.320 

7.99 

3 

10.739 

7.45 

4 

11.148 

6.96 

5 

11.557 

6.53 

6 

12.128 

5.96 

7 

12.670 

5.46 

8 

13.236 

5.00 

9 

13.927 

4.48 

10 

Run  No.  2 

15.222 

6.94 

11 

15.448 

6.14 

12 

15.659 

5.46 

13 

15.878 

4.86 

14 

16.050 

4.47 

15 

16.165 

4.15  _ 

16 

Run  No.  3 

16.685 

6.50 

16.800 

5.98 

i 

16.919 

5.53 

19  1 

17.064 

5.01 

20  '  1 

17.203 

4.54 

21 

Run  No.  4 

17.633 

6.53 

22 

17.696 

6.01 

23 

17.770 

5.50 

24 

17.849 

5.01 

25 

17.923 

4.51 

26 

Run  No.  5 

19.228 

!  7.46 

27  ! 

19.489 

7.03 

28  : 

19.789 

6.54 

29 

20.147 

!  6.01 

30  : 

20.508 

5.53 

•  31  ; 

20.858 

j  5.08 

32  ; 

21.364 

4.51 

33 

Run  No.  6 

22.545 

7.50 

34 

; 

22.784 

6.96 

35 

i 

i 

23.256 

5.97 

36 

23.487 

5.53 

37 

j - 

23.815 

4.97 
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Table  1.  Crack  Length  and  Stress  Intensity  Range  for  Several 
Fatigue  Tests  in  Specimen  No.  8353  (R  =  0.1)  (Concluded) 


File 

No. 

Test  Run 
Number 

Crack  Length 
(mm) 

Stress  Intensity 
Range,  AK,  (MPaVm) 

38 

Run  No.  7 

24.326 

9.11 

39 

25.065 

9.52 

40 

25.880 

10.00 

41 

26.695 

10.53 

42 

27.404 

11.02 

43 

28.026 

11.49 

44 

28.628 

11.97 

45 

29.246 

12.51 

46 

29.746 

12.97 

47 

30.259 

13.49 

48 

30.770 

14.03 

49 

31.173 

14.48 

50 

31.684 

15.10 

51 

31.989 

15.49 

52 

32.390 

16.03 

53 

33.109 

17.06 

54 

33.409 

17.53 

55 

i 

33.708 

18.03 

56 

34.326 

19.10 

57 

34.834 

20.08 

58 

35.243 

20.92 

59 

35.756 

22.06 

60  ! 

36.159 

23.05 

61  : 

36.568 

24.12 

62  ; 

36.871 

24.96 

63  : 

37.274 

26.18 

64  i 

37.584 

27.15 

65  1 

37.899 

|  28.22 

66 

38.105 

I  28.94 

67 

! 

38.407 

30.07 

68 

1 

! 

38.824 

31.77 

69 

39.327 

34.02 

70 

39.736 

36.04 

71 

40.472 

40.21 

72 

41.118 

44.68 
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Table  2.  Crack  Length  and  Stress  Intensity  Range  for  Several 
Fatigue  Tests  in  Specimen  No.  8354  (R  =  0.8) 


File 

No. 

Test  Run 
Number 

Crack  Length 
(mm) 

Stress  Intensity  Range, 
AK,  (MPaVm) 

1 

Run  No.  1 

10.927 

4.52 

2 

11.554 

4.01 

3 

12.266 

3.49 

4 

13.063 

2.99 

5 

14.008 

2.50 

6 

14.661 

2.20 

7 

Run  No.  2 

15.400 

4.05 

8 

15.705 

3.49 

9 

16.017 

3.01 

10 

16.391 

2.53 

11 

Run  No.  3 

17.140 

3.97 

12 

17.315 

3.53 

13 

17.562 

3.01 

14 

17.866 

2.47 

15 

Run  No.  4 

18.519 

3.45 

16 

18.908 

4.05 

17 

19.195 

4.57 

18 

19.393 

4.95 

19 

19.682 

5.58 

20 

19.870 

6.07 

21 

20.076 

6.58 

22 

Run  No.  5 

21.580 

4.47 

23 

22.075 

4.03 

24 

22.748 

3.51 

25 

1  i 

23.406 

3.05 

26  ! 

24.323 

2.52 

27  i 

! 

24.930 

2.22 

28 

Run  No.  6 

25.756 

4.98 

29  1 

26.299 

5.52 

30  ! 

26.794 

6.02 

31 

27.181 

6.47 

32 

j  27.645 

7.05 

33 

27.950 

7.46 

34 

28.329 

8.00 

35 

:  28.628 

8.45 

36 

|  29.002 

9.06 

37 

|  29.294 

9.56 

38 

i  29.583 

10.10 

39 

!  29.804 

10.50 

40 

30.099 

11.09 

41 

30.307 

11.54 

42 

30.526 

11.98 

43 

!  30.736 

12.45 

44 

!  30.963 

12.93 

45 

!  31.166 

13.58 

46 

31.440 

13.99 

47 

31.722 

14.60 
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Crack  Growth  Direction 


100  tun _  |  100  |im 


NM-312581-3 

Figure  3.  Contrast  and  gray-scale  elevation  images  of  an  area  on  the  fracture  surface  of  Specimen  8353. 
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FAST  FOURIER  TRANSFORM  (FFT)  ANALYSIS 

A  fast  Fourier  transform  (FFT)  analysis  was  applied  to  a  row  of  data  points  parallel  to  the 
crack  growth  direction,  and  the  elevation  power  spectrum  density  distribution  was  calculated. 

Each  row  contained  600  data  points;  however,  only  the  central  512  (=2y)  points  were  analyzed. 

The  EPSD  curves  for  individual  rows  exhibit  significant  differences,  as  shown  in  Figure  4.  due  to 
local  differences  in  microstructure.  Such  microstructural  effects,  although  important  from  a 
materials  science  viewpoint,  obscure  the  effect  of  loading.  Noting  that  the  stress  intensity  range, 
AK.  used  to  represent  the  loading  condition  is  an  average  value  over  the  length  of  the  crack  front 
and  not  a  local  AK  acting  on  a  single  grain  or  a  few  grains,  we  depressed  local  microstructural 
effects  and  emphasized  the  effects  of  loading  condition  by  averaging  the  EPSD  results  over  2000 
lines  wavelength-by-wavelength. 

The  averaged  EPSD  curves  from  the  different  locations  on  the  fracture  surface  are  closely 
parallel  between  wavelengths  of  0.5  pm  and  20  pm  (see  Figure  5,  for  example).  Choosing  the 
square  root  of  the  EPSD  [SQRT(EPSD)]  value  at  a  wavelength  of  5  pm  as  a  representative  index 
of  the  roughness  of  each  area  and  plotting  these  values  as  a  function  of  position  on  the  fracture 
surface,  we  obtain  the  results  shown  in  Figures  6(a)  and  6(b)  for  Specimens  8353  (R  =  0.1)  and 
8354  (R  =  0.8),  respectively.  Also  plotted  are  the  AK  values  corresponding  to  the  fracture  surface 
locations  or  crack  lengths. 

At  first  inspection,  the  behaviors  of  SQRT(EPSD)  and  AK  as  a  function  of  crack  length 
seem  quite  different,  and  a  correlation  between  load  and  roughness  is  not  clear.  For  example, 
when  the  load  ..hedding  tests  were  initiated  at  a  higher  AK,  the  SQRT(EPSD)  values  did  not  start  at 
a  higher  value,  but  instead  gradually  increased  to  a  higher  value  before  decreasing  with  decreasing 
AK.  In  the  monotonically  increasing  AK  fatigue  test,  the  SQRT(EPSD)  values  initially  decreased 
before  increasing. 

In  seeking  to  unify  these  results,  we  plotted  all  of  the  SQRT(EPSD)  values  in  Figure  6  as  a 
function  of  AK.  The  results  in  Figure  7  show  three  regions  for  each  specimen:  ( 1 )  a  linear 
relationship  in  the  low  AK  region  (less  than  about  8  MPaVm  for  8353  and  about  6  MPaVm  for 
8354).  although  a  significant  scatter  exists;  (2)  a  transition  region  (between  8  and  16  MPaVm  for 
8353  and  between  6  and  9  MPaVm  for  8354);  and  (3)  another  linear  relationship  region  at  higher 
AK  (above  16  MPa Vm  for  8353  and  9  MPa Vm  for  8354). 
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(a)  EPSD  Curve  of  Data  along  Line  400 


(b)  EPSD  Curve  of  Data  along  Line  800 
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(c)  EPSD  Curve  of  Data  along  Line  1200  (d)  EPSD  Curve  of  Data  along  Line  1600 

NM-312581-4 

Figure  4.  EPSD  Curves  for  the  roughness  profiles  along  several  lines  shown  in  Figure  3. 

(Specimen  8353) 
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Figure  5.  Averaged  EPSD  curves  from  five  locations  on  the  surface  in  Test  Run  1  for  Specimen  8353. 
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CRACK  LENGTH  (mm) 

(b)  Specimen  8354 


NM-312581-6 


Figure  6.  Square  root  of  EPSD  values  and  stress  intensity  range,  AK, 
as  a  function  of  crack  length  for  Specimens  8353  and  8354. 


The  equations  for  these  linear  relationships  are  given  below: 


£ 


8353  low  AK  region:  SQRT(EPSD)  =  4.710  +  1.067*  AK 

8353  higher  AK  region:  SQRT(EPSD)  =  4.1 16  +  0.1077*  AK 


8354  low  AK  region:  SQRT(EPSD)  =  6.776  +  0.6553*  AK 

8354  higher  AK  region:  SQRT(EPSD)  =  -2.957  +  1 . 1 56*  AK 


The  relationship  between  SQRT(EPSD)  and  AK  is  not  simple,  but  can  be  explained  in 
terms  of  crack  tip  plastic  deformation  and  crack  closure,  as  discussed  below. 


CRACK  CLOSURE  EFFECT  ON  THE  FRACTURE  SURFACE  TOPOGRAPHY 

Two  competing  processes  are  responsible  for  the  deformation  topography  of  the  fracture 
surfaces:  (1 )  the  roughening  of  the  surface  due  to  increasing  AK  or  K^  and  (2)  the  flattening  of 
the  surface  due  to  crack  closure.*  When  the  AK  values  are  low,  Figure  7  indicates  that  the 
SQRT(EPSD)  value  increases  with  increasing  AK.  In  this  region,  the  effect  of  roughening  of  the 
surface  with  increasing  AK  is  dominant. 

However,  when  AK  reaches  a  certain  value,  the  SQRT(EPSD)  starts  to  decrease  as  AK 
continues  to  increase  (transition  region  in  Figure  7).  In  this  region,  crack  closure  begins  to  have  an 
effect,  flattening  asperities  on  the  fracture  surfaces  and  gradually  overcoming  the  effect  of  surface 
roughening  with  increasing  AK.  The  surface  flattening  effect  outpaces  surface  roughening  as  AK 
increases.  Thus,  the  SQRT(EPSD)  values  decrease  despite  the  increase  in  AK. 

As  AK  continues  to  increase,  the  effect  of  surface  roughening  starts  to  overcome  the  effect 
of  surface  flattening,  and  the  SQRT(EPSD)  values  start  to  increase  again  as  AK  increases.  This  is 
the  third  region  in  Figure  7. 

With  the  above  explanation  in  mind,  we  examined  and  compared  the  results  in  Figure  6. 
Recall  that  Specimen  8353  was  tested  at  a  stress  ratio  R  =  0.1  and  8354  at  R  =  0.8,  and  that  AK 
and  K^  were  higher  for  Specimen  8353  than  for  Specimen  8354.  Since  the  AK  in  the  first  region 
is  low  and  crack  closure  is  unlikely  for  either  specimen,  the  SQRT(EPSD)  in  the  first  region  of  the 
specimen  tested  under  higher  AK  (Specimen  8353)  should  show  higher  values.  Figure  7  confirms 
this.  Conversely,  the  test  condition  with  lower  stress  ratio  should  result  in  a  higher  degree  of  crack 


‘  Crack  closure  is  the  phenomenon  of  the  faces  of  a  crack  contacting  and  perhaps  deforming  each  other  when  the  load 
is  reduced  toward  Kmn  during  load  cycling. 
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closure,  and  so  the  drop  in  the  SQRT(EPSD)  of  Specimen  8353  (with  R  =  0.1)  should  be  more 
significant  than  that  of  Specimen  8354  (with  R  =  0.8).  Figure  7  also  shows  this  to  be  the  case. 

Thus  crack  closure  can  explain  the  results  from  the  two  specimens,  so  we  examined  the 
individual  fatigue  testresults  mom  closely.  Figure  8  is  a  magnified  view  of  the  resulfr  sh°wr ,  m 
Figure  6.  In  these  load  shedding  tests,  the  AK  was  decreased  as  the  crack  extend  d.  When  AK 
above  a  certain  threshold  as  shown  in  Figure  8,  the  SQRT(EPSD)  exhibits  low  values  and  the 
higher  the  value  of  AK,  the  lower  the  SQRTfEPSD)  value.  If  the  AK  values  are  below  the 
threshold,  the  trend  in  SQRTfEPSD)  change  is  similar  to  that  of  AK. 

This  threshold  trend  is  observed  also  for  tests  in  which  AK  was  increased  as  the  crack 
grew.  Figure  8(b).  Here  the  changes  in  SQRTfEPSD)  above  the  threshold  line  are  ^  same  as  or 
the  load  shedding  test;  i.e„  the  higher  AK  is  above  the  threshold  line,  the  lower  the  SQR^ 
value.  This  trend  is  consistent  with  crack  closure  modifying  the  fracture  surfaces  at  AK  abo 

the  threshold  value. 


additional  effects  of  crack  closure 


To  delineate  the  effect  of  crack  closure  on  the  surface  roughness,  we  compared  the  EPSD 
curves  from  the  fracture  surface  areas  in  the  third  region  in  Figure  7  at  similar  AK's  of  about 
,7  MPar/m  in  die  samples  8353  (R  =  0.1)  and  8354  (R  =  0.8).  The  four  cuives  inFigure  9  (two 
for  R  =  0  1  and  the  other  two  for  R  =  0.8)  show  that,  within  the  same  R-value,  the  EPSD  curves 
have  the  same  shape  and  are  pamllel  each  other;  however,  a.  different  R-values,  dre  ^shave 
different  shapes.  The  curves  for  R  =  0.8  are  nearly  straight  between  wavelengths  of  0.5  pm  and 
20  pm;  however,  the  curves  for  R  =  0.1  deviate  from  linear  below  1  pm  and  above  4  pm.  is 
suggests  that  crack  closure  more  significandy  affects  the  features  below  the  scale  of  1 1™ 
above  4  pm.  In  dre  transition  region,  die  EPSD  values  decrease  and  the  curves  bow  more  as  AK 
increases,  as  shown  in  Figure  10.  This  suggests  that  as  AK  increases,  crack  closure  effects 

become  more  important. 

The  EPSD  changes  in  the  load  shedding  tests  were  also  examined.  We  examined  the  first 
fatigue  tes,  in  each  sample  where  crack  grow*  was  initialed  at  AK  well  above  dre  threshold I  hne 
shown  in  Figure  8.  EPSD  values  were  initially  low,  but  gradually  increased.  The  s  ape  o 
curves  during  this  change  is  shown  in  Figure  1 1 .  The  EPSD  cuives  corresponding  to  imtial  high 
AK  exhibit  bowing;  however,  as  AK  decreases  below  the  threshold  AK,  the  value  the  EPS 
curves  become  sfraighter.  Again  these  curves  suggest  that  dre  higher  AK  loading  induces  crack 

closure  effects  that  cause  bowing  of  the  EPSD  curves. 


19 


1353 


(b)  Results  for  8354 

NM-312581-7 

Possible  correlation  of  SQRT(EPSD)  behavior  with  AK  behavior. 
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Figure  9.  Comparison  of  EPSD  curves  for  tests  under  similar  stress  intensity  range  but  different  stress  ratios. 
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ELEYAnON  POWS 


We  looked  in  the  technical  literature  for  results  from  other  investigators  that  might  support 
these  findings.  Figure  12  shows  an  observation  by  Newman  of  crack  opening  stress  ratios  for 
2024-T3  aluminum  alloy  at  low  and  high  stress  ratios.14  A  striking  similarity  exists  in  the  shape  of 
the  curve  for  R  =  0.7  with  our  SQRT(EPSD)  curve  for  R  =  0.8  shown  in  Figure  8.  The  curve  in 
Figure  12  suggests  that,  due  to  the  crack  closure  and  re-deformation  of  the  fracture  topography, 
less  stress  is  necessary  to  open  the  crack.  In  other  words,  the  lower  the  opening  stress,  the  higher 
the  effect  of  surface  deformation  and  flattening  of  the  surface. 

WAVELET  ANALYSIS 

Wavelet  analysis  was  applied  to  the  topographic  data  files  by  Professor  Naoki  Saito  and 
Dr.  Jean-Marie  Aubry  of  the  Mathematics  Department  of  the  University  of  California  at  Davis 
(Appendix  A).  The  objectives  were  (1)  to  classify  the  fracture  surfaces  according  to  stress 
intensity  range,  AK,  and  (2)  to  identify  signatures  that  specify  the  loading  condition  and  can  be 

used  to  discriminate  others. 

To  conduct  the  classification  study,  the  available  data  set  needed  to  be  split  into  training  and 
test  data  sets.  The  training  data  set  was  used  to  construct  the  classification  rules,  and  the  test  data 
set  was  used  to  access  and  evaluate  whether  such  rules  can  predict  the  classes  of  the  fracture 

surfaces. 

The  Local  Discriminate  Basis15  (LDB)  algorithm  first  decomposes  the  training  signals  into  a 
basis  dictionary,  which  is  a  large  collection  of  the  basis  functions  (such  as  wavelet  packets  and 
local  cosine  functions)  localized  in  both  space  and  spatial-frequency.  Then,  it  computes 
discriminant  information  at  each  coordinate  in  the  dictionary.  The  origin?’  LDB  algorithm  uses  the 
difference  of  energy  distributions  as  the  discriminant  measure,  whereas  me  modified  LDB16  (we 
call  this  LDB+)  uses  the  probability  density  function  (pdf)  in  each  coordinate,  i.e.,  the  way  the 
coefficients  are  distributed.  Then,  a  complete  orthonormal  basis  called  LDB  (or  LDB+),  which 
“can  see”  the  distinguishing  signal  features  among  signal  classes,  is  selected  from  the  dictionary  by 
maximizing  such  discriminant  information.  After  the  basis  is  determined,  expansion  coefficients  in 
the  most  important  several  coordinates  (features)  are  fed  into  a  traditional  classifier  such  as  the 
Linear  Discriminant  Analysis  (LDA).  Finally,  the  corresponding  coefficients  of  the  test  signals  are 
computed  and  fed  to  the  classifier  to  predict  their  classes. 
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AK,  MPa-m1/2 

NM-312581-020 


Figure  12.  Crack-opening  stress  ratios 
stress  ratios.  (Ref.  14) 


for  2024-T3  aluminum  alloy  at  low  and  high 


Our  experiments  can  be  grouped  into  two  categories. 

(1)  Appnrano^o^lDl^algOTttm^on^ter^oworcolunrmoralong^a 

SSSSSSrtSr oTSKace-domain  data,  or  frequency-domam 

data). 

(2)  Application  of  2D  LDB  algorithms  to  dte  same  dam  sets  viewed  as  images. 

Summary  of  ID  Classification  Results 

.  in  general,  the  greater  the  difference  in  stress  intensity  range,  the  lower  the 
misclassification  rates. 

.  The  misclassification  rate  varies  with  the  number  of  1X>B  features  used. 

*  L^t"sffiSr“ve“e  smS"  “ 

higher  (8.4%  ~  41%). 


Summary  of  2D  Classification  Results 

.  We  computed  the  2D  power  spectrum  (modulus  of  the  2D  Fourier  tmnsform)  of 
each  image  as  the  basic  data  set. 

.  The  best  and  most  robust  results  were  obtained  by  the  LDB  based  on  the  local 
-  cosine  dictionary. 

.  Again,  the  more  different  the  stress  intensity  factors,  the  lower  the 
misclassification  rates. 

small  (20%  -  35%). 

'  SSficatomS Sght^fgtT(30%m  35°%tfferent ^  ’ 

.  The  variance  of  misclassification  rates  for  the  2D  results  is  much  smaller  than  for 
the  ID  case. 
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These  preliminary  results  suggest  that  wavelet  analysis  principles  are  well  suited  to 
overcome  the  challenges  of  extracting  and  discriminating  fracture  surface  features.  However,  the 
technique  must  be  extended  and  adapted  to  address  issues  of  noise  in  the  data  and  the  linking  of 
information  over  a  range  of  size  scales.  Moreover,  more  decisive  results  will  likely  be  obtained  by 
increasing  the  size  of  the  area  chosen  for  LDB  application.  The  results  are  discussed  in  a  later 
section,  and  suggestions  for  modifying  the  procedure  are  discussed  in  the  section  entitled  Next 

Steps. 

The  success  achieved  in  correlating  Fourier  results  with  load  parameters,  however, 
encouraged  us  to  attempt  an  analysis  of  a  field  failure.  In  the  next  section,  we  describe  how  we 
applied  the  findings  from  the  laboratory  specimens  to  interpret  the  fracture  surfaces  of  a  crack  in  an 
aircraft  turbine  component.  Appendix  B  describes  how  we  applied  fracture  surface  topographic 
analysis  to  assist  a  team  investigating  the  failure  of  a  space  shuttle  component. 


r*- 
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ANALYSIS  OF  FIELD  FAILURES 

Two  aircraft  engine  components  that  had  acquired  cracks  during  service  (Figure  13)  were 
provided  to  SRI  for  analysis.17  Component  39  had  a  part-through  crack 
Component  20  had  a  through-crack.  Ftgure  14  snows  the  racture  surfac  sjhe  " 
regions  on  both  components  am  probably  oxide  films  that  formed  during  service.  The  crack 
Component  39  apparendy  initiated  at  a  single  site  near  the  midpoint,  whereas  the  crack ' 

Component  20  o^inated  from  several  small  cracks,  which  nucleated  near  the  mtdpom  and 
merged.  The  concentric  elliptical  markings  indicate  the  crack  fronts  at  various  stages 

advance. 

topography  characterization 

we  chose  to  analyze  small  areas  of  the  crack  faces  at  regular  distances  (indicated  by  the  tick 
marks  in  Figure  15)  from  the  crack  origins  along  a  line  A-A  normal  to  the  componen  su  a  . 

opposite  component  surface.  This  area  was  divided  into  a  600  x  400  pixel  arr  y, 
of  each  pixel  was  measured.  The  data  spacing  was  0.233  pm  in  each  direction. 

Figures  16  and  17  show  a  series  of  contrast  images  and  gray-scale  topography  images. 

The  shades  of  gray  indicate  relative  height:  ligKer  areas  are  higher,  and  darker  areas  are  lower. 

The  areas  near  the  component  exteriors,  which  were  covered  by  the  film,  are  smoother  t  an 
elsewhere  and  hence  were  omitted  in  the  evaluation. 

FOURIER  ANALYSIS  AND  FRACTURE  PARAMETER  CORRELATION 

Differences  in  the  topographies  arising  from  differences  in  loading 
discerned  from  Figures  16  and  17.  To  obtain  an  alternative  representation  of  the  topography,  w 
applied  a  one-dimensional  fas,  Fourier  transfotm  operation,  row  by  row,  to  400  rews  of 
512  (ft  elevation  data  points.  Figures  18  and  19  show  the  energy  power  spectral  den  «y  (EPSD) 
curve  for  each  of  die  areas  in  Figures  16  and  17  obtained  by  averaging  die  400  rows.  In 
Figures  1 8  and  19,  the  square  mot  of  the  sum  of  die  squares  of  die  real  and  -magmaD  com^nem 
of  die  elevation  amplitudes  normalized  by  a  crack  extension  distance  ts  presented  as  function 
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Figure  13.  Two  aircraft  engine  components  that  acquired  cracks  during 
service.  (Component  39  had  a  part-through  crack  not  visible  in  this 
photo,  and  Component  20  had  a  through-crack  along  the  component 
axis.) 


NM-2037-1 


i  /  Plan  view  of  fracture  surfaces  in  Components  39  and  20. 
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Figure  15.  Fracture  surface  of  Components  39  and  20  showing  probable  crack  nucleation  sites 
(The  white  Line  A-A  and  tick  marks  indicate  the  location  of  the  areas  examined  in  detail.) 
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Figure  1 5.  Fracture  surface  of  Components  39  and  20  showing  probable  crack  nucleation  sites. 

(The  white  Line  A-A  and  tick  marks  indicate  the  location  of  the  areas  examined  in  detail.)  (concluded) 


(e)  Contrast  Image  at  570  pm  from  Surface  (f)  Elevation  Image  at  570  pm  from  Surface 


(g)  Contrast  Image  at  820  pm  from  Surface  (h)  Elevation  Image  at  820  pm  from  Surface 

NM-2037-5 


Figure  16.  A  series  of  contrast  and  gray-scale  elevation  images  at  various  locations  along 
Line  A-A  for  Component  39. 
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(i)  Contrast  Image  at  1070  pm  from 


Surface  G)  Elevation  Image  at  1 070  pm  from  Surface 


,  v  _  _ *  1C7H  i im 


from  Surface  (n)  Elevation  Image  at  1 570  pm  from  Surface 
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Elevation  Image  at  1 820  pm  from  Surface 

NM-2037-6 


Figure  16.  A  series  of  contrast  and  gray-scale  elevation  images  at  various  locations  along 
Line  A-A  for  Component  39.  (continued) 


(c)  Contrast  Image  at  320  pm  from  Surface  (d)  Elevation  Image  at  320  pm  from  Surface 


(e)  Contrast  Image  at  570  pm  from  Surface  (f)  Elevation  Image  at  570  pm  from  Surface 


(g)  Contrast  Image  at  820  pm  from  Surface  (h)  Elevation  Image  at  820  pm  from  Surface 


NM-2037-8 

Figure  17.  A  series  of  contrast  and  gray-scale  elevation  images  at  various  locations  along 
Line  A-A  for  Component  20. 
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(i)  Contrast  Image  at  1 070  pm 


(j)  Elevation  Image  at  1070  pm 


(k)  Contrast  Image  at  1 320  pm 


(I)  Elevation  Image  at  1320  pm 
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(n)  Elevation  Image  at  1570  pm 
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(o)  Contrast  Image  at  1820  pm 


(p)  Elevation  Image  at  1820  pm 

NM-2037-9 


Figure  1 7.  A  series  of  contrast  and  gray-scale  elevation  images  at  various  locations  along 
Line  A-A  for  Component  20.  (continued) 
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(q)  Contrast  Image  at  2070  pm 
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(r)  Elevation  Image  at  2070  pm 
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(s)  Contrast  Image  at  2320  pm 
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(u)  Contrast  Image  at  2570  pm 
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(t)  Elevation  Image  at  2320  pm 


(v)  Elevation  Image  at  2570  pm 
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(w)  Contrast  Image  at  2820  pm 
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(x)  Elevation  Image  at  2820  pm 
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Figure  17.  A  series  of  contrast  and  gray-scale  elevation  images  at  various  locations  along 
Line  A-A  for  Component  20.  (continued) 
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(y)  Contrast  Image  at  3050  pm  (z)  Elevation  Image  at  3050  pm 

NM-2037-1 1 

Figure  17.  A  series  of  contrast  and  gray-scale  elevation  images  at  various  locations  along 
Line  A-A  for  Component  20.  (concluded) 
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WAVE  LENGTH 
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Figure  1 8.  Elevation  powor  roectrum  density  curves  at  various  positions  along  the  Line  A-A  for  Component  39. 


WAVE  LENGTH 
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Figure  19.  Elevation  power  spectrum  density  curves  at  various  positions  along  the  Line  A-A  for  Component  20. 


spatial  frequency  or  roughness  wavelength,  tather  than  as  a  function  of  position  (as  in 
Figures  16  and  17). 

The  averaged  EPSD  curves  for  the  various  fracture  surface  areas  vtuy  in  height  and  shape 

•  u  tr^nd  with  crack  length  The  curves  are  roughly  parallel  to  one  ano 
and  show  a  characteristic  trend  with  crack  lengtn.  ffom 

“oTte^e  otigin^tiiarwere  coveredby  ,172  film,  see  Figures  20(a)  and  20(d)],  A. 
distances  k^tTlJl?  mm  from  tite  fracture  origin,  the  curves  for  the  two  components  are  nearly 

identical;  see  Figures  20b  and  20c,  for  example. 

Figure  21  shows  how  the  surface  roughness  varies  with  crack  depth.  We  arbitrarily 

ted  the  EPSD  at  a  wavelength  of  5  pm  and  plotted  the  square  root  of  this  value  versus  crac 
selected  the  because  we  observed  a  linear  relationship 

dpnth  fWe  selected  the  square  root  function  otEPbunecausewcu 

P  '  f  cpcn  value  and  AK  during  analysis  of  the  laboratory  specimens.)  If 

welnCet  ^  roughness  variation  with  cmck  lengtii  for  Component  39  as  me  variation 
of  AK  then  crack  initiation  and  early  growth  occurred  at  a  relatively  htgh  AK,  cracking  m  the 
component  interioroccurred  at  lower  AK,  and  crack  growth  closer  to  the  far  end  of  the  s^^en 
occurred  at  increasing  AK.  This  interpretation  suggests  that  the  component  failed  under  some 
modTof  bending,  because  bending  produces  stresses  that  are  highest  at  die  component  exteriors 

and  minimum  at  the  central  neutral  axis. 

■me  curve  for  Component  20  Ues  close  to  that  for  Component  39,  although  fewer  dam 
points  could  be  acquired  because  a  larger  portion  of  the  fracture  surface  w.c»ver«.^  the  oxide 
film  The  effect  of  the  oxide  film  is  to  lower  the  EPSD  values,  as  seen  in  Figure  18  and 
Figures  20(a)  and  2CM),  but  the  results  for  the  two  components  are  consent  when  no  ox.de 

exists  on  the  surface. 

Thus  the  analysis  so  far  has  provided  an  answer  to  one  aspect  of  the  load  question-the 
mode  of  loading.  The  next  question  is  the  magnitude  of  the  load.  With  die  lack  of  a  value  from 
analysis  considerations,  we  turned  to  correlation  with  darn  from  the  laboratoiy  specimens  w 

the  loading  conditions  were  known. 

Comparisons  were  made  from  two  viewpoints;  (1)  die  shapes  of  the  EPSD  curves  and 
(2)  the  heights  of  the  EPSD  curves.  The  shape  of  an  EPSD  curve  seems  to  be  sensitive  to  the 
stress-ratio,  R.  An  EPSD  curve  from  a  higher  R-value  surface  is  straighter  than  that  from  a  ower 
R-value  surface.  Figure  22  compares  EPSD  curves  of  Component  39  with  cnves  from 
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SPECTRUM  DENSITY 


laboratory  specimen  in  which  the  crack  was  propagated  at  a  stress  ratio  of  0.8.  The  shapes  of  the 
curves  match  closely;  therefore,  chose  to  correlate  the  results  of  Components  39  and  20  with  those 
from  the  specimen  tested  under  a  stress  ratio  of  0.8.  Moreover,  we  speculate  that  the  components 
experienced  a  high  R  value,  perhaps  about  0.8,  because  the  EPSD  curves  from  the  laboratory 
specimen  tested  at  R  =  0.1  were  significantly  more  convex. 

We  next  attempted  to  correlate  the  heights  of  the  curves  to  obtain  AK  values.  When  we 
correlated  the  AK  with  the  square  root  of  the  EPSD  at  a  wavelength  of  5  pm  for  the  laboratory 
specimen  tested  at  R  =  0.8,  we  obtained  the  following  relationship 

(EPSD)172  = -2.957  +  1.156  xAK  (R  =  0.8) 

Because  the  ranges  of  the  EPSD  of  both  components  were  smaller  than  those  from  the 
laboratory  specimen  (due  possibly  to  closure  effects),  we  had  to  extrapolate  to  determine  AK 
information.  The  converted  results,  shown  in  Figure  23,  suggest  that  the  cracks  in  the  engine 
components  grew  to  a  length  of  1 .6  mm  under  rather  constant  AK  conditions  of  about  5  MPa Vm 
after  initiating  at  a  higher  AK.  At  this  crack  depth,  the  stress  intensity  range  increased  strongly  and 
linearly  to  10  MPaVin  in  the  next  millimeter  of  crack  advance.  Such  a  variation  in  AK  with  crack 
depth  is  consistent  with  vibratory  bending  stresses,  because  bending  deflections  produce  higher- 
stresses  near  the  component  surface. 

WAVELET  ANALYSIS  AND  FRACTURE  PARAMETER  CORRELATION 

To  determine  whether  wavelet  analysis  could  determine  the  load  conditions  that  produced 
the  fracture  surfaces  of  Component  39,  Professor  Naoki  Saito  and  Dr.  Jean-Marie  Aubry  of  the 
Mathematics  Department  of  the  University  of  California  at  Davis  applied  the  Local  Discriminate 
Basis  algorithm  and  the  Linear  Discriminant  Analysis  algorithm  to  the  elevation  data.  Training  data 
were  obtained  from  the  laboratory  specimens  described  above,  which  had  areas  of  fracture  surfaces 
produced  at  AK  levels  of  6.6,  8.8,  1 1.0,  13.2,  and  16.5  MPaVm .  Test  data  sets  from 
Component  39  consisted  of  80  randomly  chosen  square  areas  64  pixels  on  a  side  (0.233  pm  per 
pixel).  A  majority  rule  was  then  applied  to  get  a  classification  for  the  whole  image.  When  the 
classification  rules  obtained  from  the  training  data  were  applied  to  the  test  data,  all  areas  of 
Component  39  were  indicated  to  have  been  produced  at  a  AK  of  6.6  MPaVm  or  lower.  This 
agrees  with  the  Unear  regression  results  of  Figure  23,  except  for  the  results  at  crack  lengths  larger 
than  2.0  mm.  Details  of  the  wavelet  analysis  procedures  are  given  in  Appendix  A. 
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Saito  and  Aubry  also  attempted  to  classify  the  surfaces  by  using  a  roughness  parameter 
defined  as  follows 

Rx(y)  =  JV1  +fx  (x,y)dx/Jdx 
Ry(x)  =  J  ^1 +fy2  (x,y)  dy  /  J  dy, 

where  fx  =9z/3x,  fy  =3z/3y,  and  z  =  f  (x,y)  is  a  fractograph  representing  the  elevation  data.  J dx 
represents  the  length  of  the  specimen  in  the  horizontal  direction,  which  in  fact  is  parallel  to  the 
crack  propagation  for  all  the  available  datasets.  J  dy  represents  the  length  of  the  specimen  in  the 
vertical  direction  (perpendicular  to  the  crack  propagation).  One  can  easily  see  that,  if  the  profile  is 
perfectly  flat,  then  R*  =  1  (or  Ry  =  1). 

The  variation  of  this  roughness  index  with  crack  depth  (distance  from  the  crack  origin)  for 
Component  39  is  shown  in  Figure  24.  Whether  the  parameter  is  measured  in  the  direction  of  crack 
propagation  or  parallel  to  the  crack  front,  the  trend  is  J-shaped,  in  good  agreement  with  the  trend 
for  EPSD  or  its  square  root  shown  in  Figure  21.  Furthermore,  the  roughness  parameter  was 
correlated  with  AK  by  applying  it  to  the  laboratory  specimen  and  used  to  evaluate  Component  39. 
Stress  intensity  range  values  of  1 1  MPaVrn  and  6.6  MPaVm  were  indicated  at  locations  2500 
and  1000,  respectively,  in  Figure  15(a).  These  values  agree  well  with  those  from  the  Fourier 
analysis  shown  in  Figure  23. 
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Crackfront  Location  (see  i-igi 


DISCUSSION 


This  seedling  study  has  demonstrated  a  method  for  extracting  quantitative  estimates  of  load 
parameters  flora  analysis  of  topographic  features  of  failure  surfaces.  This  is  the  first  <”*»*» 
prediction  has  been  made.  If  the  method  can  be  validated  and  made "^-friendly,  tt  wtll  represent 
a  major  advance  in  the  state  of  practice  of  failure  analysis. 

The  method  consists  of  three  parts:  (1)  characterization  and  mathematical  analysis  of 

fracture  surfaces  topography,  (2)  development  of  a  reference  datable  using  well-charactenzed 

laboratoty  specimens,  and  (3)  relating  the  reference  data  to  field  failure  topography  an 

determining  the  service  load.  The  Fourier  method  demonstrated  two  important  points:  the  shape 
of  EPSD  craves  appears  to  depend  on  the  stress-ratio,  and  the  magnitude  of  the  EPSD  (or  square 
root  of  EPSD)  relates  to  DK.  Further  exploration  of  EPSD  curve  shape  with  wavelet  techniques 
could  allow  other  fatigue  parameters  such  as  R-value  and  Kmax  to  be  determined. 


The  fast  Fourier  transfoim  analysis  of  the  data  may  not  be  optimal.  The  Fourier  results 
presented  hem  are  highly  averaged.  Therefore,  the  EPSD  curves  are  rather  smooth  and  shift  up  or 
down  within  the  0.5  to  15  nm  wavelength  range.  Individual  hues  of  topographic  dam,  however, 
show  significant  hash,  because  of  the  roughness  contribution  of  the  microstructure. 

The  topographic  relief  produced  by  the  microstructure  is  several  times  greater  than  that 
produced  by  inelastic  deformation  and  hence  masks  the  effects  produced  by  the  loads,  making  the 
extraction  of  load  parameters  difficult.  Therefore,  we  analyzed  the  averaged  EPSD  curves  to 
suppress  microstructural  effects  and  accentuate  the  load-induced  topography. 

Although  the  EPSD  curve  showed  that  features  having  certain  wavelengths  exist  on  the 
surface,  the  locations  on  the  fracture  surface  where  these  features  exist  were  not  indicated.  This  is 
a  serious  limitation  of  the  fast  Fourier  analysis  in  interpreting  fracture  surface  topograp  y. 
Therefore  we  applied  wavelet  analysis,  which  does  not  have  this  shortcoming,  to  the  elevation 
data  from  the  two  aircraft  engine  components  and  used  the  results  from  the  same  laboratory 


specimens  to  evaluate  the  results. 

To  obtain  the  training  data,  the  local  discriminant  basis  method  was  applied  to  ID  and  2D 
Fourier  treatments  of  the  elevation  data  from  the  laboratoty  specimens.  Classification  rules  were 
established  for  surfaces  produced  at  DK  values  of  6.6,  8.8, 11.0, 13.2,  and  16.5.  W  en  e  ru  es 
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were  applied  to  the  regularly  spaced  locations  on  the  fracture  surface  of  Component  39,  all  areas 
were  classified  as  DK  =  6.6.  This  result  suggests  that  nonstandard  wavelet  analysis  techniques 
will  be  required  to  discriminate  among  fracture  surfaces  produced  under  different  DK  conditions. 

We  see  three  main  reasons  wavelet  analysis  did  not  produce  more  decisive  results:  (1) 
crack  closure  effects,  (2)  the  size  of  the  area  (64  pixel  by  64  pixel)  chosen  for  LDB  analysis,  and 
(3)  noise  in  the  topography  data  file.  We  explain  how  these  aspects  influence  the  results  below. 

We  confirmed  that  crack  closure  significantly  influences  the  fracture  surface  roughness  and  leaned 
how  the  wavelet  and  FFT  results  and  their  relationship  with  DK  are  affected.  Figure  8  show  s  that 
DK  of  6.6,  8.8,  and  1 1  are  in  the  transition  regions  and  a  relationship  between  the  surface 
roughness  and  DK  in  this  region  is  not  linear.  Thus,  when  we  analyze  the  fracture  surfaces  of 
these  regions,  the  LDB  results  will  be  mixed.  However,  we  believe  that  wavelet  analysis  can 
provide  solutions. 

Crack  closure  does  influence  the  topography  of  fracture  surfaces,  and  we  believe  that  the 
degree  of  alteration  of  the  surface  is  a  function  of  the  degree  of  closure,  which  in  turn  is  a  function 
of  DK  and  the  stress  ratio.  The  crack  closure  effects  result  in  local  changes  and  thus  are  best 
characterized  with  wavelet  techniques.  We  envision  that  wavelet  can  identify  those  local  closure 
features  on  the  surface  that  can  be  used  to  classify  the  surfaces  in  the  transition  region.  One  reason 
for  the  wide  scatter  in  the  LDB  results  was  the  small  size  of  the  area  we  chose  to  investigate  relative 
to  the  material  grain  size.  With  a  pixel  spacing  of  0.233  |im,  the  area  investigated,  64  pixel  by 
64  pixel,  was  15  fim  square,  which  is  approximately  the  area  of  an  average  grain  (about 
10-20  |im).  Since  the  behavior  of  any  individual  grain  at  the  macro  crack  tip  will  differ  from  other 
grains  due  to  differences  in  orientation,  properties,  and  surrounding  grains,  the  topography  of  the 
examined  area  will  not  be  representative  of  the  material  overall.  Thus  the  local  microstructure 
dominated  the  LDB  analysis  results  and  obscured  the  effects  of  loading  condition.  It  is  necessary 
to  consider  special  way  to  assess  the  average  effect  of  loading. 

Another  source  of  error  was  noise  in  the  topography  data,  which  arose  from  the  inability  of 
the  confocal  microscope  to  get  unambiguous  light  reflection  from  high  angle  facets  in  the  fracture 
surface.  During  scanning,  the  reflected  light  intensity  was  monitored  at  every  pixel  point  as  the 
distance  between  the  fracture  surface  and  the  optical  microscope  changed.  The  light  intensity  of  a 
specific  pixel  reaches  its  peak  when  the  surface  is  in  the  focal  plane  of  the  microscope  optical 
system  and  the  microscope  records  the  value  of  the  peak  light  intensity  and  the  elevation  of  the 
pixel.  A  problem  arises,  however,  when  the  reflected  light  intensity  is  low,  because  the 
microscope  cannot  determine  when  the  peak  is  reached  at  a  specific  pixel.  Thus,  the  recorded  peak 


50 


MW.  arbitrary  and  its  elevation  is  different  from  the  surrounding  area.  This  single 
is  an  ideal  task  for  wavelet  analysis. 

Finally  the  load  shedding  experiments  used  here  are  not  die  best  type  of  experiment  for 
generating  topographie  data.  Rather,  experiments  in  which  die  stress  intensity  mcreases  as  d, 

from  the  service  failure.  Therefore,  extrapolation  was  required  to  estimate  the  service  ^ 
"e  Sanon  of  die  laboratory  reference  da.  with  the  da.  from  die  serv,«  component 

assumed  d.t  the  component ^l^^b^viOTcan  be^tssumed,  the 

DfeaUhe  differed  aacCpths  are  as  shown  in  Figrne  23.  When  superimposed  on  the  otter 
stresses^  die  component,  die  bending  stresses  should  result  in  a  profile  similar  to  that  in 

Figure  23.  .  , 

The  next  steps  are  to  develop  wavelet  and  Fourier  techniques  to  identify  and  extract  stgna  s 

from  topography  da.  that  relate  to  load  and  crack  history  parameters  and  to  build  a  da.base  of 

,  a  fnrfVkpr  cten  is  to  develop  computational  models  for  micro 

for  linking  fracture  surface  topography  features  to  load  parameters.  Future  work  .s  descnb 
the  section  entitled  Next  Steps. 
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CONCLUSIONS 


Failure-analysis-critical  load  parameter  information  is  encrypted  on  the  fracture  surfaces  of 
failed  structural  components. 

Mathematical  analysis  of  fracture  surface  topography  can  provide  quantitative  estimates  of 
the  load  conditions  that  caused  the  failure. 

A  leapfrog  advance  in  the  state  of  practice  of  failure  analysis  is  imminent. 

A  3D  fractography  analysis  technology  would  enhance  air  readiness  and  avoid  hundreds  of 
millions  in  annual  maintenance  costs. 
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NEXT  STEPS 


Much  was  learned  during  this  seedling  program.  The  load  sensitivity  of  fracture  surface 
topography  was  confirmed,  a  route  for  extracting  load  information  was  demonstrated,  the 
limitations  of  Fourier  analysis  were  shown,  xJ  the  requirements  for  extending  wavelet  analysis 
were  clarified.  We  realized  the  need  to  examine  features  at  much  higher  resolution,  the  need  for  a 
microstructural  failure  model,  the  need  to  apply  wavelet  analysis  to  larger  areas  relative  to  the 
material  grain  size,  and  the  need  to  handle  noise  in  the  elevation  data. 

required  size  scales 

Both  the  Fourier  and  the  wavelet  analyses  of  the  elevation  data  from  the  laboratory 
specimens  discriminated  fracture  surfaces  produced  under  different  load  conditions.  This  confirms 
previous  findings  suggesting  that  fracture  surface  features  are  sensitive  to  the  loading  conditions 
that  caused  the  crack.  The  physical  reason  must  he  in  the  deformation  and  microfailure  behavior  of 
the  microstructure.  Referring  to  Figure  1,  the  slip  band  activity  and  the  grain  boundary  rupture 
activity  are  two  processes  that  form  the  fracture  surfaces  in  this  nickel  alloy,  and  hence  these 
processes  must  be  sensitive  to  the  characteristics  of  the  load  history. 

The  confocal  optics  method  used  here  to  measure  fracture  surface  topography  has  a 
resolution  limit  of  0.25  pm  and  hence  is  quite  capable  of  capturing  fracture  surface  features  on  the 
scale  of  the  grains  (about  10  pm  in  the  titanium  turbine  blade  alloy,  Ti-6A1-4V).  However,  the 
steps  on  the  grain  surfaces  produced  by  interacting  slip  bands  are  typically  below  this  resolution. 

Thus,  a  higher  resolution  characterization  instrument  is  needed.  The  atomic  force 
microscope  is  capable  of  measurements  as  fine  as  several  angstroms  and  hence  capable  of 
recording  roughness  features  produced  by  slip  bands.  A  problem  then  is  combining  features  over 
a  wide  size  range.  Fortunately,  much  thought  has  been  given  to  this  problem,  and  several  schemes 
exist  for  linking  disparate  size  scales  (see  Reference  18,  for  example). 

ANALYSIS  OF  ADDITIONAL  SAMPLES  OF  DIFFERENT  STRESS  RATIOS 

In  this  program,  we  examined  the  fracture  surfaces  of  laboratory  specimens  produced 
under  two  stress  ratios,  R  =  0.1  and  R  =  0.8,  and  found  clear  differences  in  the  characteristics  of 
the  data.  To  establish  the  relationship  between  R  and  topographic  features,  we  need  to  examine 
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surfaces  produced  under  other  stress  ratios,  say,  R  =  0.3  and  R  =  0.5.  If  the  results  show  a 
consistent  trend,  then  a  unified  empirical  equation  may  predict  the  relationship  among  the 
SQRT(ESPSD),  DK,  and  stress  ratio.  Such  an  equation  can  be  used  to  extract  the  loading 
conditions  from  the  fracture  surface  topography  and  to  interpret  crack  closure  phenomena  at  the 

crack  tip. 

ROBUST  FEATURE  EXTRACTION 

From  our  classification  experiments,  we  realized  that  it  is  critically  important  to  extract  the 
features  that  are 

(1 )  Distinguishing  the  data  belonging  to  the  different  classes 

(2)  Common  and  invariant  within  the  data  belonging  to  the  same  class. 

It  is  item  (2)  that  we  need  to  address  more  for  any  tough  classification  problems  such  as  the 
fractography  data.  The  reason  for  emphasizing  the  common  features  among  the  same  class  is  both 
for  robustness  of  the  features  and  for  shedding  light  on  the  underlying  physics;  we  need  to  extract 
the  features  insensitive  to  the  variations  of  the  data  in  the  same  class. 

SPIKE  REMOVAL  BY  CONTINUOUS  WAVELETS 

The  fracture  surface  elevation  data  obtained  with  confocal  optics  microscopy  contain 
undesirable  spikes  that  may  influence  almost  all  frequency  components  in  the  Fourier  analysis. 
Since  the  high  roughness  of  the  data  contains  useful  information  that  must  be  preserved, 
conventional  denoising  techniques  cannot  be  used.  A  continuous  wavelet  transform  of  these  data 
and  the  properties  of  wavelet  coefficients  related  to  local  smoothness  of  the  data  may  be  effective  m 
detecting  and  removing  the  spikes. 

LINKING  LOCAL  FEATURES  TO  PHYSICAL  FRACTURE  MODELS 

The  specific  fracture  features  responsible  for  producing  the  "characteristic  fingerprint  of  a 
given  load  condition  in  a  given  microstructure  may  lead  to  a  much  more  realistic  fracture  model. 
That  is,  if  we  understand  how  the  microstructure  is  dissipating  energy  at  specific  load  conditions, 
the  fracture  mechanicians  and  materials  scientists  may  be  able  to  generate  a  model  that  mimics  those 
microstructural  cracking  phenomena.  This  would,  in  turn,  generate  a  model  that  links  the 
wavelet/Fourier  analysis  to  physical  reality.  One  thing  to  consider  in  this  scenario  would  be  the 
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effect  of  grain  orientation,  with  respect  to  the  crack  propagation  direction,  on  the  fractnre  features 
produced.  This  effect  could  be  used  if  or  when  we  consider  crack  propagation  within  specific, 

nonisotropic  grains. 

It  seems  to  us  that  the  newer  sets  of  dictionaries  of  bases,  such  as  local  Fourier,  brushlets, 
edgelets  and  curvelets,  which  are  capable  of  detecting  and  characterizing  oriented  and  striped 
patterns,  may  be  better  tools  to  extract  features  and  link  them  to  the  physical  models.  This  project 
requires  close  collaboration  between  mathematicians,  fracture  mechanicians,  and  matenals 

scientists. 

The  empirical  correlation  between  the  stress  intensity  and  the  EPSD  needs  to  be  understood 
in  terms  of  the  micromechanisms  that  produce  the  fracture  surface. 


Needed  at  this  point  is  a  reference  database  of  topographic  maps  and  EPSD  curves  from 
fracture  surfaces  generated  under  conditions  in  which  the  important  load  parameters  are  varied 
independently  and  systematically  over  the  entire  range  of  interest  for  each  alloy  of  interest.  Such 
databases,  known  as  fractography  atlases,  currently  exist  in  a  qualitative  form  as  a  collection  o 
SEM  micrographs  and  are  used  to  reach  subjective,  approximate  conclusions  regarding  the  causes 
of  a  failure.  An  atlas  of  topographs  and  EPSD  curves  would  allow  objective  and  quantitative 
conclusions  and  represent  a  quantum  improvement  over  current  handbooks. 
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APPENDIX  A 


CLASSIFICATION  OF  FATIGUE  FRACTURE  SURFACE  TOPOGRAPHY 
USING  MULTISCALE  AND  SPACE-FREQUENCY  METHODS 


Classification  of  fatigue  fracture  surface 
topography  using  multiscale  and 
space-frequency  methods 

Naoki  Saito  and  Jean-Marie  Aubry 
January  26,  2000 

Abstract 

We  explored  the  application  of  local  discriminant  basis  (LDB) 
methods,  which  efficiently  extracts  discriminant  local  features,  to  clas¬ 
sify  fractographs  (fracture  surface  topography  data)  of  titanium  spec¬ 
imens  according  to  the  loading  parameter  A K,  the  stress  intensity 
factor,  under  which  the  fracture  propagated. 

Our  experiments  are  mainly  grouped  into  two  categories.  1)  Appli¬ 
cation  of  ID  LDB  algorithms  along  either  row  or  column  or  circular  di¬ 
rections  with  varying  parameters  (normalizations,  filters  discriminant 
measures,  number  of  features,  space-domain  data  or  frequency- domain 
data).  2)  Application  of  2D  LDB  algorithm  to  the  same  datasets 
viewed  35  images  with  varying  parameters. 

Prom  the  results  on  the  above  experiments,  we  found  the  follow¬ 
ing.  1)  In  general,  the  larger  the  difference  between  the  A K  values  of 
these  classes,  the  lower  the  misclassification  rates.  This  suggests  that 
fracture  surface  elevation  data  definitely  contains  the  information  sen¬ 
sitive  to  A K.  2)  2D  classification  on  the  Fourier  domain  is  more  stable 
than  ID  classification,  mainly  because  we  can  assume  more  statistical 
independence  on  the  2D  image  patches,  and  the  Fourier  magnitudes 
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are  translation  invariant.  3)  It  is  still  difficult,  however,  to  single  out 
the  features  only  sensitive  to  the  change  of  A K  from  the  elevation 
data.  The  local  variability  of  the  data  is  so  severe  that  even  the  two- 
class  classification  problems  is  tough  (e.g.,  the  misclassification  rate 
ranges  around  30%  to  40%).  We  need  to  apply  the  majority  rule  on 
the  classification  results  of  many  samples  for  correct  classification. 

Based  on  these  findings,  we  recommend  the  following  actions  should 
be  taken:  1)  Develop  robust  feature  extraction  and  classification  algo¬ 
rithms,  which  are  insensitive  to  the  local  variations  in  microstructures 
2)  Link  local  features  to  physical  fracture  models  with  the  collabora¬ 
tion  with  the  material  science  experts.  3)  Build  stochastic  models  that 
can  simulate  the  similar  topographic  data  from  the  available  samples. 

4)  Integrate  multiple  sensors  and  resolutions  for  better  and  relevant 
feature  extraction  for  load  parameter  estimation. 

1  Introduction 

1.1  High  cycle  fatigue  and  structural  failure 

High  cycle  fatigue  is  a  cause  of  failure  for  materials,  especially  metals,  occur¬ 
ring  when  a  cyclic  load  is  applied  at  some  location.  L  sually  a  small  initial 
defect  already  exists,  around  which  the  stress  concentrates.  When  the  cyclic 
load  is  applied,  the  stress  amount  oscillates  between  values  Afmin  and  Kmax 
and  a  fatigue  fracture  propagates  through  the  material. 

If  this  happens  to  a  critical  part  in  an  aircraft,  such  as  a  turbine  blade, 
and  the  undetected  fractures  reach  the  point  of  final  rupture,  consequences 
can  be  dramatic.  Therefore,  a  better  understanding  of  this  phenomenon  is 
very  important.  After  the  failure  occurred,  one  would  like  to  diagnose  the 
causes  of  the  fracture  propagation,  for  instance,  the  stress  intensity  factors 
Kmin  and  Km ax  or  their  difference  A K.  This  information  is  to  be  used 
later,  for  example,  to  determine  during  which  part  of  the  flight  most  fracture 
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propagation  occurs,  and  to  improve  the  fracture-resistance  of  the  material. 

1.2  Existing  approaches 

The  method  currently  used  by  aircraft  engine  manufacturers  to  analyze  frac- 
tographs  is  by  trained  human  eyes.  These  people  usually  work  on  SEM 
images,  not  elevation  data  as  we  do. 

First  studies  at  SRI  [5,  6,  7,  8,  9,  19]  proved  some  feasibility  of  using 
the  (mean)  power  spectrum  of  adjacent  rows  of  elevation  data  as  a  tool 
to  discriminate  and  estimate  the  AJsT.  Some  aspects  of  our  work  are  an 
extension  of  this  idea  (§  5.2.3  and  5.3.2). 


1.3  Fundamental  hypotheses 

This  study  is  based  on  the  first  fundamental  assumption  that  the  history 
of  the  fatigue  parameters  is  encrypted  on  the  fracture  surface.  A  second 
hypothesis  is  that,  for  a  given  material,  only  the  A K  parameter  is  relevant 
and  governs  the  fracture  surface  topography;  all  other  variations  due,  for 
instance,  to  the  fine-scale  material  structure,  are  considered  as  random  but 
constant  in  a  probabilistic  sense.  In  other  words,  we  suppose  that  the  law  of 
the  stochastic  process  creating  the  fracture  depends  only  on  A K. 

Note  that  this  is  a  very  strong  hypothesis.  One  aspect  of  this  work  is  to 
try  to  validate  it. 


2  Mathematical  formulation 

2.1  Probabilistic  setting 

In  all  generality,  we  can  pose  the  problem  as  follows.  Suppose  that  the 
material  is  fixed,  and  let  {Q.,F,V)  be  a  probability  space.  According  to  the 
fundamental  hypotheses,  there  exists,  for  each  value  of  the  parameter  A K 
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in  a  certain  range,  a  random  process  w  €  ft  }ak{^)  whose  realizations 
represent  the  fracture  surface  topography  of  a  piece  of  material  that  was 
broke*’  under  high-cycle  fatigue  with  the  parameter  A K.  Each  realization 
is  thus  a  function  from  S  — >  R,  where  S  C  R2  is  the  area  where 
the  fracture  surface  is  analyzed.  For  a  given  location  x  €  S,  f &k{u){x)  is 
the  actual  elevation  of  this  particular  sample  at  coordinate  x  (note  that  in 
practical  applications,  5  is  a  finite  set  of  discrete  points) . 

The  problem  we  want  to  solve  is  now:  given  a  function  fx  :  S  ->  R 
representing  the  fracture  topography  of  a  sample  broken  under  some  unknown 
conditions  (X  stands  here  for  experiments),  determine  or  estimate  A K  such 
that  fx  can  be  viewed  as  a  realization  of  the  process  Jak- 

This  is  not  a  good  question  in  our  setting,  because  the  answer  may  not  be 
unique.  Rather,  we  formulate  it  in  a  Bayesian  sense:  considering  all  possible 
values  of  A K,  what  is  the  probability  density  function  (pdf)  ART  pfx  (A K) 
that  fx  Is  indeed  a  realization  of  f ak  ?  Once  this  is  answered,  we  can  chose 
“the”  A K  to  be  the  point  of  maximum  likelihood  in  pjx ,  or  its  expectation 
f  AKpfx(AK),  or  whatever  suits  best  to  the  application. 

The  calculation  of  pjx,  however,  requires  two  types  of  information:  a 
prior  probability  p0  on  A K,  and  the  law  of  the  process  u  ■-+  /**  (<*>).  The 
first  one  may  be  obtained  from  technical  considerations:  for  instance  in  the 
case  o i  an  airplane  turbine  blade,  we  may  know  what  sort  of  loads  to  expect 
in  formal  service,  and  thus  have  a  prior  probability  distribution  for  A K  (we 
can  use  uniform  distribution,  if  this  knowledge  is  not  within  the  range  of 
aeronautical  engineers’  predictions).  The  second  one  is  not  available,  as  we 
discussed  in  Introduction. 

As  is  often  the  case  in  statistics,  we  do  not  have  the  law  of  the  stochastic 
process  of  interest  and  we  have  to  estimate  it  from  a  few  sample  realizations, 
which  are  the  laboratory  samples  fractured  under  known  conditions  in  our 
case. 

Here  we  must  leave  the  world  of  the  rigorous  Bayesian  theory  for  more 
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practical  matters. 


2.2  Classification  problem 

In  our  case,  the  number  of  available  realizations  of  the  process  is  very  small 
(less  than  five  per  A K),  and  moreover  we  only  have  them  for  a  limited  set  of 
A K  values  (see  §  3).  Therefore  we  must  (temporarily)  give  up  the  estimation 
problem  and  concentrate  on  the  simpler  classification  problem. 

Let  us  give  the  classical  formulation  of  a  signal  classification  problem. 
We  first  define  appropriate  spaces  of  input  signals  (or  patterns),  extracted 
features,  outputs  (or  responses),  and  mapping  functions  among  them.  Let 
j  C  F  denote  a  signal  space  (or  a  pattern  space)  which  is  a  subset  of 
the  standard  n-dimensional  vector  space  and  which  contains  all  signals  (or 
samples/patterns)  under  consideration.  In  this  case,  the  dimensionality  of 
the  signal  space  or  equivalently  the  length  of  each  signal  is  n.  Let  2)  = 
{1,2,...  ,K}  be  a  set  of  the  class  or  category  names  to  which  the  input 
signals  belong.  We  call  this  space  a  response  space. 

Signal  classification  can  be  considered  as  a  mapping  function  (usually 
many-to-one)  d  :  3f  -*  2)  between  these  two  spaces.  Direct  manipulation 
of  signals  in  the  signal  space  for  classification  is  prohibitive  because:  1)  the 
signal  space  normally  has  very  high  dimensionality  (e.g.,  for  a  typical  segment 
of  fractograph  image  we  use  f^i  our  experiments,  n  =  64  x  64  =  4096), 
and  2)  the  existence  of  noise  or  undesired  components  (whether  random  or 
not)  in  signals  makes  classification  difficult.  On  the  other  hand,  the  signal 
space  is  overly  redundant  compared  to  the  response  space.  Therefore,  it  is 
extremely  important  to  reduce  the  dimensionality  of  the  problem,  i.e.,  extract 
only  relevant  features  for  the  problem  at  hand  and  discard  all  irrelevant 
information.  If  we  succeed  in  doing  this,  we  can  greatly  improve  classification 
performance  both  in  its  accuracy  and  efficiency. 

For  this  purpose,  we  set  up  a  feature  space  5cKk  where  k  <  n  between 
the  signal  space  and  the  response  space.  A  feature  extractor  is  defined  as  a 


map  /  :  X  -*  5,  and  a  classifier  (or  predictor )  as  a  map  g  :  5  ?)•  Let 

T  =  {(Xt>y*)}£Li  C  X  x  2)  be  a  training  (or  learning )  dataset  with  N  pairs 
of  signals  Xi  and  responses  (class  names)  y,  This  is  the  dv.taset  to  be  used  to 
construct  a  feature  extractor  /.  Let  Nc  be  the  number  of  signals  belonging 
to  class  c  so  that  we  have  N  =  Ni  +  •  •  •  +  Ata' .  Also,  let  us  denote  a  set 
of  class  c  signals  by  where  Ie  C  {1, ...  ,N}  is  a  set  of 

indices  for  class  c  signals  in  the  training  dataset  with  |/c|  =  Nc. 

Preferably,  the  performance  of  the  whole  process  should  be  measured 
by  the  misclassification  rate  using  a  test  dataset  T'  =  {(y',  (which 

has  not  been  used  to  construct  the  feature  extractors  and  classifiers)  as 
(l/N')  6(y'i  —  d(x'i)),  where  S(r  ^  0)  =  1  and  (5(0)  =  0.  If  we  use 

the  resubstitution  (or  apparent)  error  rates  (i.e.,  the  misclassification  rates 
computed  on  the  training  dataset),  we  obviously  have  overly  optimistic  fig¬ 
ures. 

For  this  application,  we  focus  on  the  feature  extractors  of  the  form 

/  =  0W  o  «P, 

where  ©(fc)  :  X  5  represents  the  selection  rule  (e.g.,  picking  most  impor¬ 
tant  k  coordinates  from  n  coordinates),  and  tf'  G  O(n),  i.e.,  an  n-dimensional 
orthogonal  matrix.  In  particular,  we  consider  matrices  representing  the  or¬ 
thonormal  bases  in  the  basis  library  (consisting  of  wavelet  packets  or  local 
trigonometric  bases)  as  candidates  for  $> .  As  a  classifier  g,  we  adopt  Linear 
Discriminant  Analysis  (LDA)  of  R.  A.  Fisher  [2]. 

2.3  Feature  extraction  using  LDB 

We  briefly  review  the  principles  of  feature  extraction  using  the  Local  Dis¬ 
criminant  Basis  algorithm  of  Coifman  and  Saito  [18,  17]. 

Let  Rn  be  the  space  of  the  signal  to  analyze.  We  assume  that  n  is  a 
dyadic  number  of  the  form  2J0,  and  we  are  given  a  hierarchical  system  of 
orthonormal  bases  for  R" . 


A-6 


Such  a  system  is  composed  of:  first  an  orthonormal  basis  B{j  which  is 
(usually)  simply  the  standard  basis  of  Rn.  Then,  Rn  is  decomposed  into  two 
orthogonal  subspaces  of  dimension  n/2,  the  first  one  having  an  orthonormal 
basis  B\  and  the  second  one  an  orthonormal  basis  B{.  The  same  decompo¬ 
sition  is  applied  to  each  of  the  subspaces,  recursively  up  to  some  maximum 
decomposition  level  ji  <  jo  (in  practice  ji  =  jo  —  1  or  jo  —  2).  At  level  j,  IT 
is  thus  decomposed  as  the  orthogonal  sum  of  2j  subspaces  of  dimension  2J0"J, 
each  of  them  with  its  orthonormal  basis  BJk  (k  =  0, . . . ,  2J  —  1).  The  number 
Nh  of  such  bases  for  R"  that  can  be  built  by  concatenating  complementary 
sub-bases  (with  maximum  decomposition  level  ji)  is  given  by  the  recurrence: 
jVj  =  Nj_ j  +  1.  These  bases  form  a  redundant  dictionary  for  R". 

For  instance,  with  ji  =  2,  the  decomposition  is  symbolized  on  Figure  1, 
and  we  have  five  possible  bases:  {Sq}.  {Bq, B}},  { Bq , ■  B%},  {B0.  Bx .  Bx }, 
{BlBlBlBl}. 


RO 

Bl 

B\ 

Bl 

B\ 

B\ 

Bl 

Figure  1:  Hierarchy  of  bases 


Such  a  hierarchical  system  of  orthonormal  bases  can  be  obtained  using  the 
discrete  wavelet  transform  (DWT)  or  local  trigonometric  transform  (Local 
Cosine,  Local  Sine,  Local  Fourier).  See  [11]  for  more  details. 

LDB  achieves  the  construction  of  a  feature  extractor  of  the  form  /  = 
@(*)  o  in  two  steps.  First,  $  is  selected  as  the  rotation  to  the  ‘‘best” 
basis  among  the  Nj}  available.  Then  0^  selects  the  k  best  vectors  in  this 
basis.  Here  “best”  is  to  be  understood  as  most  suitable  for  the  discrimination 
problem  at  hand. 

The  choice  of  the  “best”  basis  and  vectors  is  the  central  problem  of  LDB. 
For  this  we  use  some  criterion  6.  saying  that  a  basis  B  is  better  than  a  basis 
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B'  if  6{B)>6{B'). 

However,  the  number  Njy  can  be  huge,  and  it  would  be  very  inefficient  to 
check  each  possible  basis  one  by  one.  So  we  need  the  criterion  6  to  be  in  ad¬ 
ditive  form,  i.e.  for  each  vector  v  we  have  S(v)  and  6(B)  =  ^(^)*  This 

allows  the  selection  of  the  basis  to  be  done  hierarchically,  by  the  following 
procedure.  Let  B{  be  some  sub-basis  in  our  dictionary,  itself  being  decom- 
posed  in  two  sub-bases  and  If  W  >  +  then 

we  keep  {B{}  in  the  decomposition  of  Rn;  else  we  replace  it  by  the  pair  of 

bases  {B{t\B%l 

Starting  from  the  bottom  level  {B]k\k  =  0,...,2»  -  1},  this  recursion 
yields  the  basis  in  the  dictionary  that  maximizes  the  criterion  5.  Moreover, 
once  the  basis  is  selected,  the  same  criterion  (already  computed)  can  be  used 
to  select  the  k  best  vectors. 

We  have  a  few  important  remarks.  First,  imposing  additivity  is  a  huge 
restriction.  Actually,  the  true  criterion  to  use  would  be  to  minimize  the 
misclassification  rate  when  the  k  best  vectors  in  the  selected  basis  are  supplied 
to  the  classification  module  g.  But  this  criterion  is  clearly  not  additive, 
and  this  problem  creates  combinatorial  explosion.  The  choice  of  an  additive 
criterion  will  be  discussed  in  §  2.5.  Another  way  to  understand  this  remark 
is  to  remember  that  in  general,  “the  best  k  vectors  are  not  the  k  best  . 

Secondly,  once  we  resign  ourselves  to  using  an  additive  criterion,  it  is  no 
longer  necessary  to  chose  the  discriminating  vectors  in  a  basis.  We  can  simply 
take  the  k  vectors  in  the  whole  dictionary  that  maximize  6,  regardless  of  their 
linear  dependence.  However,  experimentally,  lifting  this  restriction  does  not 
improve  the  results,  but  makes  the  computation  slightly  more  demanding  (in 
memory). 

2.4  Linear  Discriminant  Analysis 

Fisher’s  LDA  first  tries  to  do  its  own  feature  extraction  by  a  linear  map 
At  :  X  5  (in  this  case  not  necessarily  orthogonal  matrix).  This  map 
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A  simultaneously  minimizes  the  scatter  of  sample  vectors  (signals)  within 
each  class  and  maximizes  the  scatter  of  mean  vectors  of  classes  around  the 
total  liean  vector.  To  be  more  precise,  let  mc:=(l/Nc)  x»  be  a  mean 
vecto-  of  class  c  signals.1  Then  the  total  mean  vector  m  can  be  defined 
as  m:=£* Li  trcmC)  where  nc  is  the  prior  probability  of  class  c  (which  can 
be  set  to  Nc/N  without  the  knowledge  on  the  true  prior  probability).  The 
scatter  of  samples  within  each  class  can  be  measured  by  the  within-class  co- 
variance  matrix  Ew:=  £*  x  ncEc,  where  Ec  is  the  sample  covariance  matrix 
of  class  c:  Ee~(l/Nc)  Eil'i (X«'C)  ~  ™c)(*|c)  -  ™c)T-  The  scatter  of  mean 
vectors  around  the  total  mean  can  be  measured  by  the  between-class  covari¬ 
ance  matrix  £(,:=£*  i  ?r c(mc  -  m)(mc  -  m)T.  Then,  LDA  maximizes  a 
class  separability  index  J  {A)~tr[{A^  EbA)~l  {AT  EWA)],  which  measures  how 
much  these  classes  are  separated  in  the  feature  space. 

This  requires  solving  the  so-called  generalized  (or  pencil-type)  eigenvalue 
problem  EbA  =  EWAA,  where  A  is  a  diagonal  matrix  containing  the  eigen¬ 
values.  Once  the  map  A  is  obtained  (normally  k  =  K  —  1),  then  the  feature 
vector  ATXi  is  computed  for  each  i,  and  finally  it  is  assigned  to  the  class 
which  has  the  mean  vector  closest  to  this  feature  vector  in  the  Euclidean 
distance  in  the  feature  space.  This  is  equivalent  to  bisecting  the  feature 
space  5  by  hyperplanes.  In  this  application  we  regard  LDA  as  a  classifier 
although,  as  explained,  it  also  includes  its  own  feature  extractor  AT.  LDA  is 
the  optimal  strategy  if  all  classes  of  signals  obey  multivariate  normal  distri¬ 
butions  with  different  mean  vectors  and  an  equal  covariance  matrix  [3],  [12]. 
In  reality,  however,  it  is  hard  to  assume  this  condition.  Moreover,  since  it 
relies  on  solving  the  eigensystem,  LDA  can  only  extract  global  features  (or 
squeezes  all  discriminant  information  into  a  few  [K  —  1]  basis  vectors)  so  that 
the  interpretation  of  the  extracted  features  becomes  difficult,  it  is  sensitive 

'The  sample  mean  operation  (1/Arc)  in  this  subsection  can  be  replaced  by  expec¬ 
tation  Ee  for  general  cases;  however,  in  this  paper,  we  focus  our  attention  on  the  cases  of 
a  finite  number  of  samples,  so  we  stay  with  the  sample  mean  operations. 
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to  outliers  and  noise,  and  it  requires  0(n3)  calculations. 


2.5  Discriminant  criterion  ■£ 

Let  us  first  concentrate  on  the  two-class  classification  problem.  Suppose  also 
that  the  probability  laws  for  each  of  these  classes  are  known:  probability  Po 
for  class  0  and  pi  for  class  1.  These  distributions  are  the  pdf’s  for  R" -valued 
random  variables  Co  and  C\ .  respectively. 

Let  v  be  a  unit  vector  in  R" .  We  consider  the  projected  probabilities  p0\v 
and  Pi|tn  that  is,  the  pdf  (on  R)  of  (Co,  v)  and  (C\ ,  v) ,  respectively. 

If  Po|t,  =  Pl |„  then  obviously  v  is  not  a  discriminating  feature:  viewed 
from  this  direction,  Co  and  Cj  look  exactly  the  same.  On  the  contrary,  if 
Po|„  is  very  different  from  pllv,  then  v  is  a  good  discriminating  feature.  Our 
criterion  8  should  measure  this  difference. 

In  information  theory',  the  difference  between  two  pdf’s  can  be  measured 
by  the  relative  entropy,  or  Kullback-Leibler  Information,  defined  by: 


J{P0\v:Pl\v)  ■=  /P0,10g 


This  formula  is  not  symmetric,  however,  so  we  may  prefer  its  symmetrized 
version: 

S(po\v,Pl\v)  :=.J(po|tMPl|r)  +  ^(Pl|f!P0|f)-  ^ 

The  symmetric  relative  entropy  can  easily  be  generalized  to  any  number 
of  classes,  taking  the  sum  of  (asymmetric)  relative  entropies  for  all  pairs  of 
classes. 

Certainly,  this  is  the  most  straightforward  choice  for  5(c).  However,  in 
practice,  we  stumble  on  the  same  ubiquitous  issue:  we  still  do  not  know  the 
probabilities  po  and  p\.  What  we  have  are  only  a  finite  number  of  realizations 
of  each  class,  constituting  the  training  set. 

Two  workarounds  are  available.  First,  one  can  try  to  evaluate  the  prob¬ 
abilities  p0|t.  and  pi|t..  Since  these  are  one-dimensional  distributions,  it  is 
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possible  to  approximate  them  by  the  histogram  of  the  values  in  the  training 
set,  or  its  smoothed  version  by  convolution  with  some  kernel  or  averaged 
with  shifted  versions  of  itself.  Averaged  Shifted  Histograms  (ASH)  are,  by 
the  way,  a  particular  case  of  kernel  convolution. 

Another  option  is  to  give  up  the  pdf-based  criteria,  and  use  the  sum  of 


energies  instead,  simply  taking 

J£(0MM  :=^l(co,«)|2log 

CO, Cl 


( l(co^)l2>\ 


where  Co  and  c\  are  taken  in  the  training  set  for  class  0  and  1,  respectively. 
And  its  symmetrized  version: 


S*(0|t;,  1\v)  :=  JB(0|t>,  l|t>)  +  Je(1|w,  0|v) 

with  the  straightforward  generalization  to  any  number  of  classes. 

As  before,  this  measures  the  difference  between  the  projection  of  each 
class  on  the  feature  v. 


2.6  Modified  Relative  Entropy 

At  some  point  in  this  study  we  tried  to  improve  the  criterion  used.  In  our 
case,  our  main  problem  is  that  we  only  have  a  very  limited  set  of  training 
data.  We  cannot  avoid  the  fact  that  they  all  come  from  the  came  laboratory 
sample,  for  instance.  As  for  the  selected  features,  we  want  them  to  be  truly 
discriminating  between  elements  of  different  A  K  classes  for  any  given  ele¬ 
vation  data,  as  long  as  the  specimen  is  the  same  metal  (e.g.  Titanium).  In 
some  sense,  we  would  like  to  reject  features  that  are  apparently  too  specific 
to  the  actual  piece  of  metal  used  for  the  training. 

Let  H  denotes  the  classical  (Shannon)  entropy  of  a  probability  distribu¬ 
tion  function.  We  recall  that  H(q)  :=  -fqlog(q)  is  small  when  q  is  a  well 
concentrated  distribution.  Suppose  that  S(po|tuPi|w)  is  I^rge  (*•£•»  very  dis¬ 
criminant).  Then  there  are  several  cases  to  which  we  need  to  pay  attention. 
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(1)  If  both  H(p 0|„)  and  H{pi\v)  are  small  (i.e.,  well  concentrated),  then 
this  vector  is  an  important  feature  extractor  for  both  classes  and  it  is  very 
unlikely  that  “chance”  (i.e.,  something  too  specific  in  this  training  data) 
would  make  this  vector  appear  important  for  both  classes.  So  we  should 
keep  this  vector. 

(2)  If,  for  instance,  H(p0\v)  is  small  and  H{p1]v)  is  large,  then  this  is 
precisely  what  we  want  to  avoid,  because  it  indicates  that  this  vector  is  good 
only  for  this  class  1  training  samples  and  may  not  carry  the  generalization 
ability  for  the  classification  problems.  LDB  might  retain  such  vectors. 

(3)  If  both  H(p0]v)  and  H(px\v)  are  large,  then  this  vector  is  specific  to 
neither  class  0  nor  class  1.  However,  there  is  no  reason  to  discard  this  vector 
since  it  still  carries  good  discriminant  information.  So,  we  keep  this  one  too. 

From  the  above  consideration,  we  propose  a  Modified  Relative  Entropy 
based  on  pdfs: 

M{po\v,Px\v)  =  S(p0\v,Pi\v)  -Amax(H{polv),H(pllv))  +  Bmm{H{p0\v),H(pl]v)) 

where  A  and  B  are  positive  constants  balancing  the  relative  importance  one 
may  want  to  impose  on  effects.  We  then  search  the  basis  that  maximizes  the 
sum  of  M  (or  its  obvious  modification  in  the  case  of  more  than  two  classes) 
from  the  dictionaries  of  bases. 

3  Data  sets 

The  following  elevation  matrices  were  provided  by  SRI.  We  mention  the  ma¬ 
terial  (a  Titanium  alloy),  a  reference  identification,  the  A K  (when  known), 
the  dimension  in  pixels,  the  horizontal  and  vertical  scale  of  the  originals 
(matrices  used  for  classification  were  rescaled  at  the  same  scale). 
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3.1  Titanium  sample  #8354 

This  is  a  laboratory  sample  that  was  submitted  to  a  stress  test  with  constant 
ratio  R  =  Kmix/Km\n  =  0.8.  Parameter  A K  —  Kmix  —  Amjn  varies  from  2.2 
to  14.6  (units:  ksix/iiT)- 

The  stress  test  is  composed  of  a  series  of  runs.  The  first  run  starts  with 
a  value  of  A K  high  enough  for  the  fracture,  initiated  by  an  initial  dent,  to 
propagate.  The  value  of  A K  is  then  gradually  reduced,  until  it  is  under  the 
threshold  below  which  the  fracture  does  not  propagate  anymore.  Another 
run  is  then  started.  The  last  run  receives  increasing  values  of  A K  until  the 
final  failure. 

3.1.1  First  data 

This  batch  of  data  is  the  first  dataset  we  received.  The  complete  runs  of 
8353  and  8354  came  later. 


Name 

A  K 

Dimension 

xy  scale 

2  scale 

(/im/pixel) 

(^m/grav  level) 

8354/Ti6.0 

6 

718  x  478 

0.36 

78.25/256 

8354/Ti6_l 

6 

479  x  638 

0.29 

62.66/256 

8354/Ti6_2 

6 

479  x  638 

0.29 

62.66/256 

8354/Ti6_3 

6 

479  x  638 

0.29 

62.66/256 

8354/Ti8_0 

8 

718  x  478 

0.36 

78.25/256 

8354/Ti8.1  j 

8 

62.66/256 

8354/Ti8_2 

8 

479  x  638 

0.29 

62.66/256 

8354/Ti8_3 

8 

479  x  638 

0.29 

62.66/256 

8354/TilO.O 

10.1 

718  x  478 

0.36 

78.25/256 

8354/TilO.l 

479  x  638 

0.29 

80.42/256 

8354/TilO_2 

1BW 

479  x  638 

0.29 

80.42/256 

8354/TilOA 

10.1 

479  x  638 

0.29 

80.42/256 

8354/Til2_0 

12 

718  x  478 

0.36 

116.59/256 

8354/Til2_l 

12 

479  x  638 

0.29 

88.11/256 

8354/Til2_2 

12 

479  x  638 

0.29 

105.80/256 

8354/Til2_3 

12 

479  x  638 

0.29 

137.08/256 

8354/Til5_0 

15 

■■ 

142.37/256 

8354/Til5_l 

15 

479  x  638 

0.29 

128.60/256 

8354/Til5_2 

15 

479  x  p38 

0.29 

84.09/256 

8354/Til5_3 

15 

479  x  638 

0.29 

128.60/256 
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3.1.2  Run  1 


Name 

A  K 

Dimension 

8354/Til0927 

4.52 

600  x  2000 

0.23 

43.01/256 

600  x  2000 

0.23 

53.95/256 

8354/Til2266 

3.49 

600  x  2000 

0.23 

58.42/256 

8354/Til3063 

2.99 

600  x  2000 

0.23 

58.42/256 

8354/Til4008 

2.50 

600  x  2000 

0.23 

58.42/256 

8354/Til4661 

0.23 

58.42/256 

3.1.3  Run  2 


Name 

A  K 

Dimension 

xy  scale 
(pm/pixel) 

z  scale 

(/um/gray  level) 

8354/Til5400 

4.05 

600  x  2000 

0.23 

65.07/256 

8354/Til5705 

0.23 

65.07/256 

8354/T16017i 

3.01 

600  x  2000 

0.23 

65.07/256 

8354/Til6391 

0.23 

72.39/256 

3.1.4  Run  3 


Name 

A  K 

Dimension 

xy  scale 
(^m/pixel) 

z  scale 

(pm/grav  level) 

8354/Til7140 

3.97 

600  x  2000 

0.23 

72.39/256 

8354/Til7315 

3.53 

600  x  2000 

0.23 

72.39/256 

8354/Til7562 

0.23 

72.39/256 

8354/Til7866 

2.47 

600  x  2000 

0.23 

72.39/256 
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3.1.5  Run  5 


Name 

A  K 

Dimension 

2  scale 

(/xm/gray  level) 

8354/Til8519 

0.23 

72.39/256 

8354/Ti 18908 

0.23 

72.39/256 

8354/Til9195 

0.23 

72.39/256 

8354/Til9393 

4.95 

600  x  2000 

0.23 

72.39/256 

8354/Til9682^ 

0.23 

72.39/256 

8354/Til9870 

0.23 

72.39/256 

8354/Ti20076 

6.58 

600  x  2000 

0.23 

72.39/256 

3.1.6  Run  7 


Name 

AK 

Dimension 

xy  scale 
(jim/pixel) 

2  scale 

(^im/grav  level) 

8354/Ti21580 

4.47 

600  x  2000 

0.23 

72.39/256 

8354/Ti22075 

4.03 

600  x  2000 

0.23 

8354/Ti22748 

3.51 

600  x  2000 

0.23 

72.39/256 

8354/Ti23406 

3.05 

600  x  2000 

72.39/256 

8354/Ti24323 

2.52 

600  x  2000 

8354/Ti24930 

0.23 

72.39/256 
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3-1.7  Run  8 


Name 

A  K 

Dimension 

xy  scale 
(/xm/pixel) 

2  scale 

(/xm/grav  level) 

8354/Ti25756 

4.98 

600  x  2000 

0.23 

72.39/256 

8354/Ti26299 

5.52 

600  x  2000 

0.23 

72.39/256 

8354/Ti26794 

6.02 

600  x  2000 

0.23 

72.39/256 

8354/Ti27181 

6.47 

600  x  2000 

0.23 

72.39/256 

8354/Ti27645 

7.05 

600  x  2000 

0.23 

84.00/256 

8354/Ti27950 

7.46 

600  x  2000 

0.23 

96.38/256 

8354/Ti28329 

8.00 

600  x  2000 

0.23 

96.38/256 

8354/Ti28628 

8.45 

600  x  2000 

0.23 

109.07/256 

8354/Ti29002 

9.06 

600  x  2000 

0.23 

109.07/256 

8354/Ti29294 

9.56 

600  x  2000 

0.23 

109.07/256 

8354/Ti29583 

10.10 

600  x  2000 

0.23 

112.71/256 

8354/Ti29804 

10.50 

600  x  2000 

0.23 

112.71/256 

8354/Ti30099 

11.09 

600  x  2000 

0.23 

112.71/256 

8354/Ti30307 

11.54 

600  x  2000 

0.23 

112.71/256 

8354/Ti30526 

11.98 

600  x  2000 

0.23 

139.67/256 

8354/Ti30736 

12.45 

600  x  2000 

0.23 

139.67/256 

8354/Ti30963 

12.93 

600  x  2000 

0.23 

139.67/256 

8354/Ti31166 

13.58 

600  x  2000 

0.23 

139.67/256 

8354/Ti31440 

13.99 

600  x  2000 

0.23 

151.97/256 

8354/Ti31722 

14.60 

600  x  2000 

0.23 

162.30/256 

3.2  Titanium  sample  #8353 

This  is  a  laboratory  sample  that  was  submitted  to  a  stress  test  with  constant 
ratio  R  —  •Kmax/^min  =  0.1.  Parameter  A K  —  Kmax  varies  from 

4.15  to  44.68. 
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3.2.1  Run  1 


Name 

A  K 

Dimension 

xy  scale 
(^m/pixel) 

z  scale 

(jzm/grav  level) 

8353/Ti09924 

8.52 

600  x  2000 

0.23 

74.93/256 

8353/Til0320 

7.99 

600  x  2000 

0.23 

74.93/256 

8353/Til0739 

600  x  2000 

0.23 

93.60/256 

8353/Tilll48 

600  x  2000 

0.23 

93.60/256 

8353/Til 1557 

6.53 

600  x  2000 

0.23 

120.22/256 

8353/Til2128l 

5.96 

600  x  2000 

0.23 

120.22/256 

8353/Til2670 

5.46 

600  x  2000 

0.23 

120.22/256 

8353/Til3236 

600  x  2000 

0.23 

120.22/256 

8353/Til3927 

4.48 

600  x  2000 

0.23 

120.22/256 

3.2.2  Run  2 


Name 

A  K 

Dimension 

xy  scale 
(//m/pixel) 

z  scale 

(^rn/grav  level) 

8353/Til5222 

0.23 

120.22/256 

8353/Til5448 

6.14 

600  x  2000 

0.23 

8353/Til5659 

5.46 

600  x  2000 

0.23 

120.22/256 

600  x  2000 

0.23 

120.22/256 

835?/Til6050 

4.47 

600  x  2000 

0.23 

120.22/256 

8353/Til6165 

4.15 

600  x  2000 

0.23 

120.22/256 
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3.2.1  Run  1 


Name 

A  K 

Dimension 

xy  scale 
(^m/pixel) 

z  scale 

(jxm/gray  level) 

8353/Ti09924 

8.52 

600  x  2000 

0.23 

74.93/256 

8353/Til0320 

0.23 

74.93/256 

8353/Til0739 

7.45 

600  x  2000 

0.23 

93.60/256 

8353/Til 1148 

6.96 

600  x  2000 

0.23 

93.60/256 

8353/Til 1557 

6.53 

600  x  2000 

0.23 

120.22/256 

8353/Til2128 

5.96 

600  x  2000 

0.23 

120.22/256 

8353/Til2670 

5.46 

600  x  2000 

0.23 

120.22/256 

8353/Til3236 

5.00 

600  x  2000 

0.23 

120.22/256 

8353/Til3927 

4.48 

600  x  2000 

.  0.23 

120.22/256 

3.2.2  Run  2 


Name 

A  K 

Dimension 

xy  scale 
(jxm/pixel) 

z  scale 

(/xm/gray  level) 

8353/Til5222 

6.94 

600  x  2000 

0.23 

120.22/256 

600  x  2000 

0.23 

120.22/256 

8353/Til5659 

5.46 

600  x  2000 

0.23 

120.22/256 

8353/Til5878 

4.86 

600  x  2000 

0.23 

120.22/256 

mi 

600  x  2000 

0.23 

120.22/256 

8353/Til6165 

0.23 

120.22/256 
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3.2.3  Run  3 


Name 

A  K 

Dimension 

xy  scale 
(//m/pixel) 

2  scale 

(/xm/gray  level) 

8353/Til6685 

6.50 

600  x  2000 

0.23 

120.22/256 

8353/Til6800 

5.98 

600  x  2000 

0.23 

120.22/256 

8353/Til6919 

5.53 

600  x  2000 

0.23 

120.22/256 

8353/Til7064 

5.01 

600  x  2000 

0.23 

120.22/256 

8353/Til7203 

4.54 

0.23 

120.22/256 

3.2.4  Run  4 


Name 

A  K 

Dimension 

xy  scale 
(//m/pixel) 

z  scale 

(/xm/gray  level) 

8353/Til7633 

0.23 

120.22/256 

8353/Til7696 

6.01 

0.23 

120.22/256 

600  x  2000 

0.23 

120.22/256 

8353/Til7849 

5.01 

600  x  2000 

0.23 

120.22/256 

600  x  2000 

0.23 

120.22/256 

3.2.5  Run  6 


Name 

A  K 

Dimension 

xy  scale 
(/zm/pixel) 

z  scale 

(/xm/gray  level) 

600  x  2000 

0.23 

120.22/256 

600  x  2000 

0.23 

120.22/256 

8353/Til9789 

6.54 

600  x  2000 

120.22/256 

8353/Ti20147 

6.01 

600  x  2000 

■ 

120.22/256 

8353/Ti20508 

5.53 

600  x  2000 

8353/Ti20858 

5.08 

600  x  2000 

0.23 

120.22/256 

8353/Ti21364 

0.23 

120.22/256 
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F 


3.2.6  Run  7 


Name 

A  K 

Dimension 

xy  scale 
(/im/pixel) 

z  scale 

(ym / gray  level) 

8353/Ti22545 

7.50 

600  x  2000 

0.23 

120.22/256 

8353/Ti22784 

6.96 

0.23 

120.22/256 

8353/Ti23256 

5.94 

600  x  2000 

120.22/256 

8353/Ti23487 

5.53 

600  x  2000 

■a 

120.22/256 

8353/Ti23815 

4.97 

600  x  2000 

0.23 

120.22/256 

3.2.7  Run  8 


Name 

A  K 

Dimension 

xy  scale 
(/xm/pixel) 

z  scale 

(^m/gray  level) 

8353/Ti24326 

9.11 

600  x  2000 

0.23 

120.22/256 

8353/Ti25065 

9.52 

600  x  2000 

0.23 

120.22/256 

8353/Ti25880 

10.00 

600  x  2000 

0.23 

120.22/256 

8353/Ti26695 

10.53 

600  x  2000 

0.23 

120.22/256 

8353/Ti27404 

11.02 

600  x  2000 

0.23 

120.22/256 

8353/Ti28026 

11.49 

600  x  2000 

0.23 

120.22/256 

8353/Ti28628 

11.97 

600  x  2000 

120.22/256 

8353/Ti29246 

12.51 

600  x  2000 

120.22/256 

8353/Ti29746 

12.97 

600  x  2000 

0.23 

120.22/256 

8353/Ti30259 

IHBM 

600  x  2000 

0.23 

120.22/256 

8353/Ti30770 

600  x  2000 

0.23 

120.22/256 

8353/Ti31173 

14.48 

600  x  2000 

0.23 

120.22/256 

8353/Ti31684 

0.23 

120.22/256 

8353/Ti31989 

15.49 

600  x  2000 

0.23 

120.22/256 

8353/Ti32390 

0.23 

120.22/256 
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3.3  Component  #39  (titanium) 

Component  #39  is  one  of  the  aircraft  engine  components  that  actually  failed 
in  service.  The  A K  parameter  is  unknown  in  this  case.  Fractographs  were 
taken  at  nine  regularly  spaced  locations,  starting  from  the  edge  of  the  com¬ 
ponent  where  the  crack  initiated. 


Name 

A  K 

Dimension 

xy  scale 
(/zm/ pixel) 

z  scale 

(/xm/gray  level) 

B39/TiEDG 

? 

719  x  478 

0.23 

29.32/256 

B39/Ti0250 

? 

719  x  478 

0.23 

29.32/256 

B39/Ti0500 

? 

719  x  478 

0.23 

29.32/256 

B39/Ti0750 

? 

719  x  478 

0.23 

29.32/256 

B39/Til000 

? 

719  x  478 

0.23 

29.32/256 

B39/Til500 

?  1 

719  x  478 

0.23 

29.32/256 

B39/Til750 

? 

719  x  478 

0.23 

29.32/256 

B39/Ti2000 

? 

719  x  478 

0.23 

29.32/256 

B39/Ti2500 

? 

719  x  478 

0.23 

29.32/256 

4  Data  preparation 

As  one  can  see,  it  is  very  difficult  to  conduct  statistical  analysis  on  these 
datasets,  for  we  have  a  very  small  number  of  samples  (maximum  five  images 
in  each  A K  class)  of  high  dimension  (up  to  600  x  2000).  It  would  make  no 
sense  to  try  to  extract  features  in  this  space.  Instead,  we  simulate  a  large 
number  of  samples  in  a  smaller  space  by  sub-sampling  the  images.  This 
artificially  creates  a  nice  statistical  set,  but  one  must  keep  in  mind  that  it 
is  basically  biased  by  the  fact  that  the  samples  are  not  at  all  independent: 
they  come  from  the  same  piece  of  metal,  sometimes  they  even  overlap.  One 
must  be  very  careful  on  this  point,  and  try  to  correct  the  bias  as  much  as 
possible. 
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Also,  one  must  check  that  the  horizontal  and  vertical  scales  of  the  files 
we  want  to  compare  are  the  same,  and  rescale  them  if  needed. 

4.1  Cne-dimensional  sampling 

We  can  sample  rows  or  columns  of  the  matrix  (the  “row”  direction  corre¬ 
sponding  to  the  direction  of  propagation  of  the  fracture),  or  along  other 
shapes  such  as  circles.  Optionally,  we  could  take  their  Fourier  transform  and 
study  them  in  the  frequency  domain. 

One-dimensional  computation  is  easier  and  faster,  and  it  was  the  first  way 
we  tried  to  validate  the  method.  However,  we  should  keep  a  critical  eye  on 
the  results,  because  the  dependence  between  the  samples  in  this  case  is  very 
strong  (by  continuity,  two  neighboring  rows  have  about  the  same  profile). 

Moreover,  at  the  very  beginning  of  this  study,  we  were  given  only  five 
images,  each  one  in  its  A K  class  (8354/Ti6.0,  8354/Ti8.0,  8354/Til0_0, 
8354/Til2_0,  8354/Til5.0).  So  we  had  to  split  them  to  make  the  training 
and  test  data  sets:  first  half  of  the  matrix  in  the  training  set,  the  second  one 
in  the  test  set,  eventually  with  some  randomization.  In  this  case,  the  test  set 
is  very  dependent  from  the  training  set,  and  the  significance  of  the  results 
may  be  very  weak.  For  instance,  the  classification  could  be  (and  actually 
was)  based  on  the  mean  elevation  or  slope  of  the  sample,  which  depend  only 
on  the  d^La  acquisition  procedure  and  not  on  the  A K.  This  is  not  what  we 
want:  we  want  to  select  features  that  are  pertinent  only  to  the  parameter 
A K.  We  tried  to  correct  this  baleful  effect  by  normalizing  the  mean  and  the 
first  order  moment  to  zero,  so  that  each  sample  has  the  same  overall  elevation 
and  slope. 

Later,  other  locations  (but  still  from  the  same  piece  of  metal)  for  each 
class  were  added  to  our  collection  (8354/Ti AA_7t,  ti  =  1,2,3).  Since,  in  each 
class,  locations  1  and  2  do  not  overlap,  we  could  use  one  of  them  as  training 
set,  the  other  one  as  test  set,  for  a  more  realistic  trial. 
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4.2  Two-dimensional  sampling 

We  also  prepared  and  examined  two-dimensional  samples.  We  took  small 
square  “patches”  of  size  32  x  32  or  64  x  64  (powers  of  two  are  more  suitable  to 
multiresolution  analysis).  Usually  they  were  taken  at  random,  but  sometimes 
with  a  fixed  pattern  such  as  checker-board  tiling.  Optionally,  we  could  take 
their  Fourier  transform  and  study  them  in  the  frequency  domain.  A  typical 
number  of  patches  is  128  on  each  image. 

In  the  case  of  two-dimensional  sampling,  the  problem  of  the  dependence 
between  neighboring  samples  seems  less  present.  We  do  not  have  the  means 
to  give  a  precise  formulation  of  this  statement,  because  this  requires  a  proba¬ 
bilistic  model  of  the  surface,  which  we  do  not  have,  but  a  coarse  justification 
is  given  by  the  remark  that,  when  using  a  regular  tiling,  a  square  of  side  c 
only  has  4c  points  that  are  close  to  points  of  its  neighbors,  while  its  total 
number  of  points  c2. 

4.3  Computational  aspects 

The  computations  were  done  on  an  Intel  Pentium  II  machine  (450  MHz  CPU, 
128  MB  RAM)  running  Linux  as  an  operating  system. 

The  data  preparation  and  feature  extraction  routines  were  implemented 
jn  C++.  The  files  were  converted  and  manipulated  as  Scientific  Data  Sets  in 
Hierarchical  Data  Format  (a  file  *ormat  developed  by  the  National  Center  for 
Supercomputing  Applications,  freely  available  at  http://hdf.ncsa.uiuc.edu). 

For  one-dimensional  analysis,  the  Local  Discriminant  Basis  code  (ldb) 
and  its  version  using  Averaged  Shifted  Histograms  (ldbkash)  had  already 
been  developed  by  one  of  us  (Saito). 

For  two-dimensional  analysis,  the  program  to  compute  the  energy  dis¬ 
tribution  on  the  basis  dictionary  (tfenergy2d),  and  the  program  using  the 
output  of  the  latter  to  compute  the  LDB  (ldb2d)  had  also  already  been 
written.  The  code  computing  the  Averaged  Shifted  Histograms  of  the  coef- 
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ficients  in  the  basis  dictionary  (tfkash2d),  as  well  as  the  code  using  these 
histograms  to  compute  the  probability  distributions-based  LDB  (ldbkash2d) 
were  developed  as  a  part  of  this  project. 

The  actual  classification  using  LDA  was  done  in  S-plus  using  a  code  from 

Venables  and  Ripley  [20]. 

5  Overview  of  the  results 

5.1  Note  on  the  interpretation 

The  results  are  given  in  the  form  of  misclassification  rates,  that  is,  the  per¬ 
centage  of  total  errors  made  by  the  classifier  over  the  total  number  of  sam¬ 
ples.  The  classifier  was  first  tested  on  the  training  data  (resubstitution  error 
rate),  then  on  the  test  data.  We  often  give  both  in  the  following  tables,  but 
as  explained  in  Introduction,  the  second  number  is  the  one  that  has  more 
significance. 

For  two-class  classification  problems,  the  worst  misclassification  rate  should 
be  50%,  in  which  case  the  classifier  operates  “at  random” ,  and  is  of  no  use. 
Misclassification  rates  >  50%,  if  they  occur,  indicate  a  potential  flaw  in  our 
fundamental  hypotheses. 

Misclassification  rates  <  50%  are  positive  results.  Though  an  error  rate 
of,  say,  35%,  may  not  seem  very  encouraging,  we  should  remember  that  we 
are  trving  to  classify  the  whole  images,  not  only  the  small  samples  that  were 
created  for  doing  the  statistics. 

So  we  can  apply  a  rule  deciding  that  the  class  of  the  whole  image  is  the 
class  of  the  majority  of  the  samples.  Assuming  that  we  have  N  independent 
sample  on  the  image,  and  that  the  probability  of  making  an  error  on  one 
sample  is  p,  the  probability  pst  of  having  a  majority  of  errors  is 

m  fN\ 

*»-  E  (*)  a -**>"-*. 
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Using  the  Central  Limit  Theorem,  we  can  approximate  this  probability 
for  large  N  by 


where  erfc  is  the  normal  tail  error  function  defined  by 

erfc(l):=7^i  e~t/2dt 

As  long  as  p  <  1/2,  the  error  rate  for  the  whole  image  decreases  expo¬ 
nentially  with  N. 

5.2  One-dimensional  classification 

5.2.1  Data  preparation:  one  file  per  class 

As  we  saw  it  before,  constructing  a  statistical  classifier  requires  a  large  num¬ 
ber  of  training  samples.  In  our  case,  the  number  of  available  images  (el¬ 
evation  matrices)  is  very  limited;  actually,  when  we  started  this  stud y,  we 
only  had  one  image  for  each  of  five  values  of  A K:  8354/Ti6_0,  8354/Ti8.0, 
8354/TilO.O,  8354/Til2.0,  and  8354/Til5.0. 

Since  a  preliminary  study  had  shown  some  results  using  one-dimensional 
Fourier  transform  on  the  rows  of  these  images,  we  began  our  exploration  in 
the  same  direction.  For  the  simplest  problem  of  two-class  classification  (for 
instance  between  8354/Ti6.0  and  8354/Til5-0),  we  would  take  half  of  the 
rows  of  each  image  for  training,  and  the  second  half  for  testing  the  classifier. 

Other  options  include  :  taking  columns,  randomized  rows  or  columns, 
random  circles...  We  also  tried  to  take  off  the  mean  and/or  the  first  order 
moment  of  each  sample,  in  order  to  suppress  the  artifact  caused  by  the  overall 
elevation  of  the  samples.  Here  is  an  example  of  misclassification  rates  between 
8354/Ti6.0  and  8354/Til5.0,  obtained  using  ldb  with  Daubechies  10  filter  to 
extract  the  50  most  discriminant  (in  the  sense  of  symmetric  relative  entropy) 
features. 


A-26 


The  first  number  represents  the  misclassification  rate  when  the  classifier 
is  applied  to  the  training  data  itself;  the  second  number  is  for  the  test  data. 


raw 

mean  0 

first  moment  0 

both  0 

rows 

0%  30% 

0%  1% 

0%  31% 

1.5%  23% 

columns 

0%  0% 

0%  0% 

0%  1.8% 

0%  1.1% 

random  rows 

0%  0% 

0%  0% 

0.42%  0.84% 

2.1%  3.3% 

random  col. 

0%  0% 

0%  0% 

0.14%  0% 

0.42%  0.7% 

random  circles 

3.1%  6.1% 

0.63%  4.2% 

12%  14% 

35%  44% 

Unfortunately,  these  figures  are  not  statistically  significant,  for  the  simple 
reason  that  training  and  test  samples  are  not  statistically  independent:  they 
come  from  the  same  image!  This  bias  is  even  worse  when  they  are  random¬ 
ized.  Moreover,  we  see  that  when  the  mean  and  the  first  order  moment  are 
removed,  the  performance  drops.  This  suggests  that  classification  is  actually 
done  on  the  overall  elevation  and  slope  of  each  image,  rather  than  on  its 
intrinsic  surface  characteristics. 

Nevertheless,  this  kind  of  very  primitive  experiment  can  yield  some  inter¬ 
esting  conclusion.  We  explored  further  the  two-class  one-dimensional  classi¬ 
fication  on  all  possible  pairs  of  the  five  classes  we  had  at  that  point.  Here 
are  the  results  obtained  using  ldb  with  Daubechies  10  filter,  50  most  dis¬ 
criminant  features,  random  columns,  mean  and  first  order  moment  set  to 
zero.  Results  above  the  diagonal  use  the  first  half  of  randomized  columns  for 
training,  second  half  for  testing;  results  below  the  diagonal  use  the  second 
half  of  randomized  columns  for  training,  first  half  for  testing. 


AK  =  6 

AK  =  8 

AK  =  10 

AK  =  12 

AK  =  15 

AK  =  6 

8.8%  13% 

AK  =  8 

6.8%  11% 

12%  17% 

2.8%  5.7% 

0.7%  0.42% 

A  K  =  10 

4.3%  5.8% 

8.8%  12% 

4.2%  6.4% 

3.2%  3.8% 

AK  =  12 

7.2%  7.2% 

4%  4.5% 

4.5%  6.3% 

0.14%  0.42% 

ART  =  15 

0.14%  0.84% 

0.56%  0.97% 

1.8%  2.5% 

0%  0.56%  1 

A-27 


As  previously  noticed,  because  of  the  statistical  dependence,  these  results 
are  completely  biased  (in  regard  of  what  would  be  the  performance  of  the 
classifip-  in  a  ceal  situation).  But  a  general  tendency  can  be  observed:  the 
more  different  the  A Ks,  the  lower  the  misclassification  rate.  This  is  an 
important  confirmation  of  our  fundamental  hypothesis. 

5.2.2  Two  files  per  class 

A  more  realistic  test  of  our  methods  is  of  course  the  case  when  the  samples 
used  to  test  the  classifier  are  independent  from  the  ones  used  to  train  it. 
In  order  to  keep  all  other  parameters  (the  physical  conditions  of  the  experi¬ 
ments)  as  neutral  as  possible  for  the  classification  problem,  training  and  test 
samples  were  taken  from  the  same  piece  of  metal  at  different  locations  (A K 
being  the  same  within  each  class).  These  samples  are  noted  8354/TiA K.n, 
with  Atf=6,  8,  10,  12,  15  and  n=l,  2,  3  corresponding  to  three  different 
locations.  Locations  1  and  2  do  not  overlap  at  all,  and  thus  we  can  use  them 
for  the  test  where  training  and  test  samples  are  independent. 

Here  we  show  the  results  for  Ti6/Til5  classification,  with  various  sam¬ 
pling  schemes,  using  Daubechies  10  wavelets  for  ldb  or  ldbkash  and  extracting 
50  basis  vectors.  Ml  means  that  the  mean  and  the  first  moment  are  set  to 
zero.  Here  the  first  number  is  the  misclassification  rate  obtained  with  files 
#1  as  craining  set,  files  #2  as  test  set;  the  second  number  is  obtained  with 
files  #2  as  training  set,  files  #1  as  test  set. 


ldb 

Ml  +  ldb 

ldbkash 

rows 

100%  50% 

30%  50% 

100%  50% 

columns 

9.5%  41% 

8.4%  40% 

19%  37% 

rows,  randomly  shifted 

8.1%  50% 

48%  47% 

7.7%  42% 

columns,  randomly  shifted 

23%  43% 

49%  46% 

24%  42% 

random  circles 

75%  50% 

47%  45% 

96%  50% 

Bad  results  and  extreme  variation  suggests  that  ldb(kash)  is  much  too  sen¬ 
sitive  to  the  particular  feature  of  the  samples.  The  asymmetry  between  each 
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pair  of  results  is  quite  surprising  :  in  almost  all  cases  the  pair  Ti6.2/Til5.2 
seems  better  for  training  than  Ti6.1/Til5.1.  This  could  result  from  some 
special  property  of  these  files,  b”  ,  this  is  obviously  not  what  we  are  looking 

for. 

Figure  3  shows  the  10  most  discriminant  LDB  vectors  computed  for  the 
Ti6.1/Til5.1  classification  problem.  Some  of  the  vectors  suggests  that  the 
discriminant  information  may  be  in  the  Fourier  (i.e.,  wavenumber  or  spatial 
frequency)  domain.  This  observation  lead  us  to  the  next  set  of  experiments. 


200  *00  600  800  1000 


Figure  3:  Example  of  top  10  seleted  basis  vectors 


5.2.3  Fourier  domain 

We  then  conducted  the  classification  experiments  in  the  Fourier  domain,  and 
use  a  basis  of  Local  Cosines  to  analyze  it.  This  is  a  refinement  of  the  original 
idea  of  looking  at  the  power  spectrum  of  each  sample:  the  effect  of  using  the 
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LCT-based  LDB  is  to  (automatically)  find  a  discriminating  portion  of  the 
Fourier  domain. 

We  FFT-ed  each  row,  column  and  circle,  and  comput'd  the  energy  dis¬ 
tributions  in  the  Fourier  domain.  These  data  were  supplied  to  the  LDB 
algorithm.  Note  that  when  no  filter  was  employed,  which  means  that  LDA 
classification  was  directly  applied  on  the  whole  energy  distributions,  the  com¬ 
putation  time  severely  increased  because  we  used  all  (or  half)  of  the  512 
coordinates. 


dlO+s 

let 

Ida  all 

Ida  2nd  half 

rows 

24%  41% 

18%  35% 

20%  45% 

51%  51% 

columns 

14%  39% 

10%  41% 

5.9%  36% 

0%  50% 

random  circles 

16%  50% 

67%  50% 

65%  48% 

0%  50% 

The  last  column  shows  the  result  obtained  when  we  tried  the  classifi¬ 
cation  based  on  the  second  half  of  the  Fourier  domain  (high  frequencies). 
This  attempt  was  suggested  by  other  evidence  from  SRI  that  some  useful 
information  may  exist  in  the  highest  frequency  region. 

However,  the  results  are  still  not  very  conclusive  (slightly  better  than  in 
the  spatial  domain,  but  not  much). 

5.2.4  Summary  of  results  for  one-dimensional  classification 

•  In  general,  the  more  different  the  stress  intensity  factors  (A K),  the 
lower  the  misclassification  rates. 

•  The  misclassification  rate  varies  with  the  number  of  LDB  features  used 
(violating  independence  of  each  sample). 

•  If  both  the  training  and  test  datasets  are  selected  from  the  same  region 
randomly,  then  the  misclassification  rates  are  very  small  (0.4  %  to  17 

%). 
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•  If  the  training  and  test  datasets  are  selected  from  the  two  different  re¬ 
gions  without  randomization,  then  the-  misclassification  rates  are  much 
higher  (8.4  %  to  41  %)- 

5.3  Two-dimensional  classification 

We  thus  decided  to  move  on  to  two-dimensional  classification.  Although 
theoretically  the  problems  are  very  similar,  in  practice  a  number  of  technical 
difficulties  appear.  The  main  one  is  that  the  dimensionality  of  the  signal 
space.  In  the  one-dimensional  problems,  the  dimension  of  the  signal  space 
was  the  length  of  the  rows  or  columns,  at  most  700  points.  In  two-dimensional 
problems,  if  we  want  to  take  patches  of  reasonable  size  as  our  samples,  say 
64  x  64,  the  dimension  jumps  to  4096.  Here  the  need  for  extracting  the 
significant  features  before  doing  the  classification  is  even  bigger,  because  the 
classifiers  such  as  LDA  perform  very  badly  (and  need  a  huge  computation 
time)  in  high  dimensions. 

We  work  here  on  the  same  files  as  before:  8354/TiAA"_n,  with  A K 
6,8,10,12,15  and  n  =  1,2,3  corresponding  to  three  different  locations.  In 
the  following  we  used  location  1  as  training  samples  and  location  2  as  test 
samples,  and  vice-versa. 

On  each  image,  we  take  a  certain  number  (128)  of  small  square  samples 
(64  x  64)  as  our  training  or  test  datasets.  Here  the  patches  are  randomly 
sampled.  Note  that  in  these  two-dimensional  experiments,  we  did  not  remove 
the  mean  and  the  first  moment  from  each  sample. 

5.3.1  Spatial  domain 

The  following  results  were  obtained  using  ldb2d  with  Daubechies  10  filter, 
extracting  the  50  most  discriminant  features.  First  number  of  each  pair  is 
the  training  misclassification  rate,  the  second  one  is  the  test  misclassification 
rate.  Results  above  the  diagonal  use  file  #1  for  training,  file  #2  for  testing; 
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Results  below  the  diagonal  use  file  #2  for  training,  file  #1  for  testing. 


A  K  =  6 

AK  =  8 

AK  =  10 

AK  — 12 

AK  —  15 

AK  =  6 

18%  29% 

21%  45% 

21%  44% 

6%  46% 

AK  =  8 

23%  30% 

29%  59% 

23%  55% 

7%  44% 

AK  =  10 

19%  47% 

26%  57% 

35%  52% 

8%  57% 

AK=  12 

17%  46% 

24%  55% 

28%  54% 

6%  63% 

AK  =  15 

7%  48% 

8%  39% 

6%  55% 

6%  64% 

5.3.2  Fourier  domain 

We  obtained  the  misclassification  rates  using  the  50  most  discriminant  fea¬ 
tures  computed  by  ldb2d  with  LCT  on  the  energy  distributions  of  the  image 
patches  in  the  Fourier  domain.  First  number  of  each  pair  is  the  training  mis¬ 
classification  rate,  the  second  one  is  the  test  misclassification  rate.  Results 
above  the  diagonal  use  file  #1  for  training,  file  #2  for  testing;  Results  below 
the  diagonal  use  file  #2  for  training,  file  #1  for  testing. 


AK  =  6 

AK  =  8 

AK  =  10 

AK  =  12 

AK  =  15 

AK  =  6 

16%  48% 

15%  48% 

8%  29% 

2%  36% 

AK  =  8 

27%  49% 

16%  43% 

12%  30% 

4%  31% 

AK  =  10 

17%  50% 

18%  39% 

AK  =  12 

11%  34% 

12%  28% 

23%  47% 

13%  52% 

AK  =  15 

2%  33% 

2%  30% 

The  results  appear  to  be  slightly  better  in  the  Fourier  domain.  Figure  4 
shows  the  typical  Fourier  domain  energy  distributions  on  these  two  classes 
of  surfaces,  and  the  partition  patterns  of  the  Fourier  domain  computed  by 
the  LDB  with  LCT. 

5.4  Two-dimensional  classification  using  new  data 

After  doing  the  experiments  described  in  §  5.3,  we  received  the  complete 
elevation  data  for  samples  8353  and  8354.  Among  the  different  runs,  A K 
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A K  =  6  (Fourier)  A K  =  15  (Fourier)  Localization  (LCT) 

Figure  4:  Fourier-transformed  samples  (energy)  and  the  selected  LCT  local¬ 
ization 


takes  a  large  range  of  values,  that  can  be  used  later  for  the  regression  problem 
of  estimating  A K  from  samples.  However,  we  still  concentrate  on  the  simpler 
two-class  classification  problem  here,  the  classes  being  chosen  among  A K  = 
6, 8, 10, 12, 15,  for  the  purpose  of  comparison  with  the  previous  experiments. 
In  ea'\.  run,  we  take  the  closest  values  to  these  AK’s: 

•  in  the  A K  =  6  class  we  have:  8354/Til9870  and  8354/Ti26794  for  sam¬ 
ple  8354,  and  8353/Til2128, 8353/Til6800, 8353/Til7696,  8353/Ti2014 
8353/Ti23256  for  sample  8353; 

•  in  the  A K  =  8  class  we  have:  8354/Ti28329  for  sample  8354,  and 
8353/Til0320  for  sample  8353; 

•  in  the  A K  =  10  class  we  .have:  8354/Ti29583  for  sample  8354,  and 
8353/Ti25880  for  sample  8353; 
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•  in  the  A K  =  12  class  we  have:  8354/Ti30526  for  sample  8354,  and 
8353/Ti26628  for  sample  8353; 

•  in  the  A K  =  15  class  we  have:  8354/Ti31722  for  sample  8354,  and 
8353/Ti31684  for  sample  8353. 

This  time  we  can  use  the  two  different  pieces  as  training  data  and  test 
data,  which  are  now  completely  independent.  For  different  results,  we  use 
here  ldbkash2d  with  the  modified  relative  entropy  presented  in  §  2.6.  The 
test  is  run  in  the  Fourier  domain  using  LCT,  and  the  number  of  features 
selected  is  5.  Results  above  the  diagonal  use  8353  for  training,  8354  for 
testing;  Results  below  the  diagonal  use  8354  for  training,  8353  for  testing. 


AK  =  6 

AK  =  8 

AK  =  10 

AK  =  12 

AK  =  15 

AK  =  6 

16%  33% 

16%  33% 

16%  33% 

15%  41% 

AK  =  8 

33%  18% 

44%  44% 

43%  56% 

42%  50% 

AK  =  10 

33%  16% 

40%  53% 

47%  54% 

45%  53% 

AK  =  12 

32%  19% 

41%  54% 

44%  47% 

45%  49% 

AK  =  15 

32%  20% 

30%  58% 

42%  53% 

47%  49% 

With  a  few  exceptions,  the  results  are  clearly  negative  (the  same  ex¬ 
periments  with  the  classical  relative  entropy  are  neither  better,  nor  worse). 
Several  reasons  can  be  imagined,  but  we  certainly  need  to  review  our  fun¬ 
damental  hypotheses.  Most  likely,  the  parameter  A K  is  not  the  only  one  to 
influence  the  surface  topography  in  this  case,  and  we  must  take  into  account 
other  factors,  such  as  the  ratio  R  =  Kmax/Kmin,  which  is  different  on  the 
samples  8353  (R  -  0.1)  and  8354  (R  =  0.8).  Some  other  evidence  from  SRI, 
using  the  mean  Fourier  power  spectrum  of  rows  in  these  samples,  corroborate 
this  hypothesis. 

5.5  Classification  of  Component  #39  data 

We  finally  conducted  a  series  of  experiments  on  the  fractographs  taken  on 
Component  #39  (see  page  3.3),  one  of  the  aircraft  engine  components  that 
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actually  failed  in  service.  The  A K  parameter  in  this  case  is  truly  unknown. 
Since  we  did  not  develop  any  regression  code,  the  best  thing  we  could  do  was 
to  apply  to  each  fractograph  the  classifier  constructed,  03  explained  above, 
to  discriminate  between  two  known  A K  values. 

Training  samples  were  taken  from  the  8353  data  set  in  a  first  experiment, 
and  from  both  8353  and  8354  data  sets  in  a  second  experiment.  We  used  all 
files  corresponding  to  the  values  of  A K  previously  analyzed;  6,  8,  10,  12,  15. 
The  same  parameters  as  in  §5.4  were  used  for  the  LDB/LDA  algorithm.  As 
the  test  datasets,  80  samples  of  size  64  x  64  were  randomly  chosen  in  each 
C39  fractographs,  and  a  majority  rule  was  then  applied  in  order  to  get  a 
classification  for  the  whole  image. 

Using  the  classification  rules  for  Ti6/TiX  where  X=8,10,12,15,  all  the  C39 
fractographs  were  classified  as  Ti6,  with  an  overwhelming  majority.  In  all 
other  classification  experiments,  such  as  Ti8/Til0,  the  classification  results 
appeared  random,  and  each  case  the  majority  vote  was  close  to  50  %. 

The  results  seem  to  indicate  that  the  C39  fractographs  were  fractured 
under  the  stress  conditions  lower  than  DK  <  6.  This  observation  is  also 
supported  by  the  following  experiments  using  the  roughness  parameters.  Let 
z  —  f(x,  y)  be  a  fractograph  representing  the  elevation  data.  We  quantified 
the  roughness  of  a  profile  of  a  fractograph  by 

Rx{y)  =  j  v/l  +  /J(x,j/)dx/  J dx 
Ry(x)  =  J  ^/l  +  /92(x,y)  dy  /  J d y, 

where  fx  =  dz/dx,  fy  =  dz/dy.  J  dx  represents  the  length  of  the  specimen 
in  the  horizontal  direction,  which  in  fact  is  parallel  to  the  crack  propagation 
for  all  the  available  datasets.  /  dy  represents  the  length  of  the  specimen 
in  the  vertical  direction  (perpendicular  to  the  crack  propagation).  One  can 
easily  see  that  if  the  profile  is  perfectly  flat,  then  Rx  =  1  (or  Ry  =  1). 

Figure  5  compares  the  roughness  parameters  of  the  specimens  of  8354 
Run  8  (and  one  specimen  from  Run  5)  with  those  of  C39  specimens. 
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Figure  5:  Comparison  of  roughness  parameters  of  profiles  of  fractographs. 
(a)  8354  Run  8  Rx ;  (b)  8354  Run  8  Ry\  (c)  C39  Rx;  (d)  C39  Ry.  For  (a) 
and  (b),  the  numbers  in  the  vertical  axis  represent  DK  of  the  specimen.  For 
(c)  and  (d),  the  numbers  in  the  vertical  axis  represent  the  location  of  the 
specimen  in  /xm.  The  box  in  each  plot  represents  the  50  %  of  the  population 
of  the  roughness  parameters  for  the  given  fractograph.  The  dot  in  the  center 
of  the  box  represents  the  median  of  the  roughness  parameters  for  the  given 
specimen.  The  circles  represent  extremely  high  values  such  as  outliers. 
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6  Conclusion 

We  summarize  our  study  on  these  fractography  datasets  as  follows: 

•  For  the  two-class  classification  problems,  in  general,  the  larger  the  dif¬ 
ference  between  the  A K  values  of  these  classes,  the  lower  the  mis- 
classification  rates.  This  suggests  that  fracture  surface  elevation  data 
definitely  contains  the  information  sensitive  to  A K. 

•  2D  classification  on  the  Fourier  domain  is  more  stable  than  ID  classi¬ 
fication,  mainly  because  we  can  assume  more  statistical  independence 
on  the  2D  image  patches,  and  the  Fourier  magnitudes  are  translation 
invariant. 

•  It  is  still  difficult,  however,  to  single  out  the  features  only  sensitive  to 
the  change  of  A K  from  the  elevation  data.  The  local  variability  of  the 
data  is  so  severe  that  even  the  two-class  classification  problems  is  tough 
(e.g.,  the  misclassification  rate  ranges  around  30%  to  40%).  We  need 
to  apply  the  majority  rule  on  the  classification  results  of  many  samples 

'  for  correct  classification  as  described  in  §  5. 

7  Recommendations 

Based  on  our  experience  on  the  above  experiments,  we  recommend  that  the 

following  actions  should  be  taken. 

1.  Develop  robust  feature  extraction  and  classification  algorithms. 

2.  Develop  an  algorithm  to  remove  undesired  spikes. 

3.  Link  local  features  to  physical  fracture  models. 

4.  Build  stochastic  models. 

5.  Integrate  multiple  sensors  and  resolutions. 
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We  now  describe  each  project,  which  can  be  further  divided  into  several 
subprojects. 

7.1  Robust  Feature  Extraction  and  Classification 

From  our  classification  experiments,  we  realized  that  it  is  of  critical  impor¬ 
tance  to  extract  the  features  that  are:  1)  Distinguishing  the  data  belonging 
to  the  different  classes,  and  2)  Common  and  invariant  within  the  data  be¬ 
longing  to  the  same  class.  It  is  2)  that  we  need  to  address  more  for  any 
tough  classification  problems  such  as  the  SRI  fractography  data.  The  reason 
for  emphasizing  the  common  features  among  the  same  class  is  both  for  ro¬ 
bustness  of  the  features  and  for  a  possibility  of  shedding  a  light  to  a  better 
understanding  of  the  underlying  physics:  we  need  to  extract  the  features 
insensitive  to  the  variations  of  the  data  in  the  same  class. 

Both  the  original  and  the  probabilistic  versions  of  the  local  discriminant 
basis  (let  us  call  them  LDB  and  LDB+,  respectively)  methods  tried  to  em¬ 
phasize  more  about  the  aspect  1)  above,  although  the  criteria  used  in  LDB+ 
implicitly  and  weakly  used  the  aspect  2).  In  our  projects,  we  will  more 
explicitly  explore  the  aspect  2). 

7.1.1  Development  of  LDB++ 

W~  are  currently  implementing  and  refining  the  discriminant  measure  we 
described  in  §  2.6.  Let  H(p0\v),  H(pi\v)  be  the  entropy  of  class  0  and  class  1 
samples  projected  on  a  given  basis  vector  v  in  the  dictionary  of  bases  such 
as  the  local  Fourier  dictionary  or  wavelet  packet  dictionary.  Let  S(po|t»Pi|<;) 
be  the  symmetric  relative  entropy  that  is  used  for  the  LDB-t-  method.  From 
the  discussion  in  §  2.6,  we  propose  to  use  the  following  new  criteria: 

M(po\v,Px\v)  =  ^(poi^PiiJ  -  >lmax(H(po|r),/f(pi|v))  +  Bm\n{H(p0\v),  H(pllv)). 
Another  possibility  is: 

A/ (po|v! Pi|t>)  =  5(p0|f,Pi|t>)  -  A\H{p0\v)  -  H{pi\v)\ ■  {H{po\v)  +  H{pi\v)). 
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A  and  B  are  positive  constants  balancing  the  relative  importance  one  may 
want  to  impose  on  effects.  We  then  search  the  basis  that  maximizes  the 
sum  of  M(p0\v,pi\v)  from  the  dictionaries  of  bases.  In  the  experiment  we 
conducted  in  §  5.4,  we  only  tested  the  one  of  the  above  formula  with  A  = 
B  =  0.5.  We  need  to  conduct  more  tests  and  research  how  to  find  optimal 
A  and  B. 

7.1.2  Subtraction  of  common  structures  using  the  joint  best  basis 

We  should  also  discard  the  truly  common  and  shared  coordinates  between  the 
two  classes  for  any  classification  tasks  since  they  are  useless  for  that  task.  To 
do  this  more  explicitly,  we  have  reached  a  simple  idea  of  common  structure 
subtraction  using  the  joint  best  basis  algorithm.  This  strategy,  although  very 
simple,  must  be  tried  for  any  practical  classification  problems. 

Step  0:  Expand  all  the  signals  of  both  classes  into  a  dictionary  of  bases. 

Step  1:  Find  a  joint  best  basis  that  is  efficient  to  compress  both  classes. 

Step  2:  Starting  from  the  most  energetic  coordinates,  mark  the  coordinates 
that  give  us  low  discriminant  values  measured  by  relative  entropy. 

Step  3:  Remove  such  energetic  but  non-discriminant  coordinates  from  both 
classes  of  signals. 

Step  4:  Apply  the  LDB,  LDB+,  or  LDB-I— I-  for  the  residuals. 

This  strategy  should  really  enhance  the  genuine  difference  between  these  two 
classes. 

7.1.3  The  local  Fourier  dictionary,  local  shifts,  and  image  segmen¬ 
tation 

Another  major  development  we  would  like  to  incorporate  is  the  local  Fourier 
dictionary  as  our  basic  dictionary  for  feature  extraction.  The  local  Fourier 
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dictionary  has  many  attractive  properties  that  the  wavelet  packets  and  local 
cosine  dictionaries  cannot  have:  the  local  shift  invariance,  the  oriented  fea¬ 
ture  extraction  capability,  to  name  a  few  (see  e.g.,  [16]).  These  properties  are 
particularly  attractive  for  the  fracture  surface  classification  because:  1)  the 
corresponding  features  within  the  same  regions  should  be  aligned  properly 
before  the  LDB  feature  extraction;  2)  we  expect  that  the  surface  images  with 
larger  magnifications  will  reveal  the  striation  patterns  locally  depending  on 
the  stress  conditions.  To  treat  the  feature  alignment  problems,  we  recently 
proposed  the  idea  of  the  Best  Local  Translation  Basis  (BLTB)  using  the  local 
Fourier  dictionary  [16].  At  this  point,  this  algorithm  is  designed  for  ID  sig¬ 
nals.  We  will  implement  and  improve  this  algorithm  for  2D  images.  For  the 
oriented  feature  extractions  and  more  generally,  image  segmentation  based 
on  the  oriented  patterns,  we  recently  launched  investigation  of  segmentation 
via  sparsity  constraints  in  the  local  Fourier  dictionary.  The  idea  here  is  to 
fully  utilize  the  best  basis  algorithm  using  the  sparsity  constraint  to  capture 
the  homogeneously  and  isotropically  oriented  regions  as  much  as  possible. 
We  expect  that  the  regions  containing  more  than  two  different  orientations 
should  be  split  into  finer  boxes.  Therefore,  the  idea  is  to  detect  the  smallest 
boxes  in  the  best  basis  partition  patterns  to  identify  the  texture  boundary, 
and  the  largest  boxes  as  a  part  of  the  homogeneous  texture  regions,  and  then 
merge  the  smaller  boxes  into  the  larger  boxes  if  they  share  the  same  texture. 
These  algorithms,  i.e.,  the  BLTB  and  the  texture  segmentation  via  sparsity 
constraints  can  be  significant  by-product  that  may  be  generated  from  this 
research  since  it  is  easy  to  foresee  such  algorithms  can  be  useful  for  the  other 
context  that  is  of  DARPA’s  interest  (e.g.,  ATR,  other  pattern  recognition 
applications). 

7.1.4  Robust  entropy  estimation 

We  have  also  realized  that  the  robust  entropy  estimation  from  a  given  set  of 
data  is  critically  important  for  all  of  the  projects  above.  Currently,  we  are 


using  the  histogram-based  entropy  estimator  of  Hall  and  Morton  [4]  and  its 
variants.  These  entropy  estimators  are  good  if:  1)  the  number  of  available 
samples  are  large;  2)  the  true  pdf  satisfies  a  certain  regularity  conditions 
(e.g.,  no  heavy  tails,  etc,);  and  3)  the  underlying  distribution  is  ID.  There¬ 
fore,  if  one  of  these  conditions  are  violated,  we  need  to  consider  alternatives. 
The  local  Fourier  dictionary  immediately  violates  the  condition  3)  since  their 
expansion  coefficients  are  complex  numbers  (i.e.,  2D  in  nature).  Therefore, 
we  need  good  and  robust  entropy  estimator  for  the  local  Fourier  dictionary. 
Even  for  the  local  cosine  and  wavelet  packet  dictionaries,  the  condition  1) 
and  2)  may  be  violated  from  time  to  time.  The  entropy  estimators  we  are 
currently  examining  are  kernel-based  density  estimator  and  graph-theoretical 
density  estimator  (e.g.,  the  one  based  on  the  minimum  spanning  trees).  With 
our  Ph.D.  student,  Jen-Jen  Lin,  we  recently  obtained  a  method  to  compute 
entropy  without  using  the  pdf  [10].  This  method  uses  the  higher  order  cumu- 
lants  instead  of  the  pdf’s.  This  also  works  for  the  multidimensional  distribu¬ 
tions.  We  have  a  plan  to  examine  its  stability  and  compare  its  performance 
with  those  of  the  pdf-based  methods. 

7.2  Spike  Removal  by  Continuous  Wavelets 

As  we  already  mentioned,  the  fractography  elevation  data  contains  the  un¬ 
desirable  spikes  that  may  influence  almost  all  frequency  components  in  the 
Fourier  analysis.  Since  the  roughness  of  the  data  contains  useful  information, 
which  must  be  preserved,  conventional  denoising  techniques  such  as  wavelet 
thresholding  of  Donoho  and  Johnstone  cannot  be  employed;  those  cannot 
distinguish  rough  features  from  spikes.  We  suggest  the  use  of  a  continuous 
wavelet  transform  of  these  data,  and  the  properties  of  wavelet  coefficients 
related  to  local  smoothness  (the  so-called  Holder  regularity)  of  the  data,  to 
detect  and  remove  the  spikes. 
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7.3  Linking  Local  Features  to  Physical  Fracture  Mod¬ 
els 

This  is  one  of  the  key  component  of  our  future  project,  and  will  be  con¬ 
ducted  with  the  fractography  and  material  science  experts.  If  we  understand 
the  specific  fracture  features  responsible  for  producing  the  “characteristic 
fingerprint”  of  a  given  load  condition  in  a  given  microstructure,  this  may 
lead  to  a  much  more  realistic  fracture  model.  That  is,  if  we  understand  how 
the  microstructure  is  dissipating  energy  at  specific  load  conditions,  the  frac¬ 
ture  mechanists  and  material  scientists  may  be  able  to  generate  a  model  that 
mimics  those  microstructural  cracking  phenomena.  This  would,  in  turn,  gen¬ 
erate  a  model  that  links  the  wavelet/Fourier  analysis  to  physical  reality.  One 
thing  to  consider  in  this  scenario  would  be  the  effect  of  grain  orientation, 
with  respect  to  the  crack  propagation  direction,  on  the  produced  fracture 
features.  This  could  be  used  if/when  we  consider  crack  propagation  within 
specific,  non-isotropic,  grains. 

It  seems  to  us  that  the  newer  sets  of  dictionaries  of  bases,  such  as  local 
Fourier,  brushlets  [13],  and  edgelets  [1],  which  are  capable  of  detecting  and 
characterizing  oriented  and  striped  patterns,  may  be  ideal  tools  to  extract 
features  and  link  them  to  the  physical  models.  This  project  requires  real 
close  collaboration  with  fracture  mechanists  and  material  scientists. 

7.4  Building  Stochastic  Models 

While  the  fracture  mechanists  and  material  scientists  develop  the  relevant 
fracture  models,  we  also  would  like  to  build  stochastic  models  that  are  able 
to  generate  as  many  fracture  surfaces  as  one  wishes,  which  are  similar  to  the 
observed  fracture  surfaces  of  a  given  loading  condition.  If  we  can  successfully 
build  such  models  for  different  loading  conditions,  then  the  model  parameters 
should  reflect  the  difference  of  the  loading  conditions.  This  in  turn  may  be 
incorporated  to  the  physical  fracture  models  of  Section  2.3.  We  have  a  couple 
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of  projects  in  this  directions. 

7.4.1  Simple  stochastic  models  using  the  least  statistically-dependent 
bases 

Statistical  independence  is  one  of  the  most  desirable  properties  as  a  coordi¬ 
nate  system  in  terms  of  redundancy  reduction.  If  the  coordinates  are  sta¬ 
tistically  independent,  modifying  one  coefficient  does  not  affect  the  other 
coefficients  at  all.  Therefore,  if  we  can  find  such  a  coordinate  system  for  a 
given  class  of  images,  we  can  use  this  system  for  a  variety  of  purposes,  e.g., 
compression,  modeling,  and  simulation.  In  reality,  however,  it  is  difficult  to 
obtain  such  a  coordinate  system  because:  1)  a  given  image  may  not  con¬ 
sist  of  independent  components  in  the  first  place;  or  2)  even  so,  numerical 
algorithms  to  obtain  such  a  basis  become  extremely  complicated  and  expen¬ 
sive.  Based  on  these  observations,  we  very  recently  developed  an  algorithm 
that  quickly  selects  a  basis  (from  the  dictionaries)  that  is  “closest”  to  the 
statisticallv-independent  coordinate  system  in  the  sense  of  relative  entropy 
[14],  [15].  The  simplest  stochastic  model  of  fracture  surfaces  is  to  assume 
that  the  LSDB  coordinates  are  statistically  independent,  and  each  fracture 
surface  obeys  the  stochastic  process  that  can  be  specified  by  the  product 
of  empirical  marginal  pdf’s  of  the  LSDB  coordinates.  We  can  then  easily 
simulate  and  sample  new  images  from  this  model  by:  1)  independently  sam¬ 
pling  the  LSDB  coefficients  using  the  standard  sampling  methods  such  as 
the  inversion  method  or  the  rejection  method;  and  2)  performing  the  inverse 
transform  from  the  LSDB  coordinates  to  the  pixel  coordinates. 

7.4.2  Stochastic  models  incorporating  dependence  among  local 
features 

If  the  LSDB  coordinates  do  not  provide  us  with  truly  or  nearly  independent 
coordinates,  we  will  incorporate  the  dependence  among  the  LSDB  coordi¬ 
nates.  W7e  will  start  with  a  simple  dependence  model  such  as  the  pairwise 
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dependence  model,  then  develop  more  elaborate  models  such  as  the  Markov 
chains  or  Markov  Random  Fields  on  the  local  features/LSDB  coordinates. 

The  simplest  dependence  model  is  the  pairwise  dependence  model  that 
approximates  the  true  pdf  by  a  product  of  bivariate  pdf  s.  The  issue  of 
this  simple  model  then  is  how  to  make  pairs  of  the  LSDB  coordinates.  We 
will  use  the  mutual  information  to  quantify  the  dependence  between  the  two 
coordinates.  Once  we  make  the  pairs  of  the  coordinates,  we  can  sample 
and  simulate  the  coefficients  using  the  standard  simulation  techniques  such 
as  the  rejection  method.  By  looking  at  the  simulation  results,  we  know 
how  important  the  dependence  structure  is  for  the  texture  models.  For  the 
coordinates  of  complex  numbers,  such  as  the  local  Fourier  coordinates,  we 
can  also  compute  the  mutual  information  on  the  magnitude  and  the  phase 
separately,  which  is  extremely  interesting  to  examine. 

If  the  pairwise  dependence  is  not  enough  to  model  the  images  faith¬ 
fully,  we  will  proceed  to  build  more  elaborate  probability  models  such  as 
the  Markov  chain  on  the  LSDB  coordinates.  Here,  the  simplest  algorithm  is 
to  start  from  the  most  energetic  coordinate,  then  find  the  most  dependent 
coordinate  on  the  first  coordinate,  then  find  the  most  dependent  coordinate 
on  the  second  coordinate,  and  repeat  until  we  exhaust  the  indices.  If,  in  the 
middle  of  the  process,  the  dependence  of  all  the  remaining  coordinates  on 
the  already-selected  coordinates  becomes  too  weak,  then  we  can  break  the 
chain  at  this  point,  and  th^n  start  the  new  chain  from  the  most  energetic  co¬ 
ordinate  that  is  still  left.  In  both  the  simple  pairwise  model  and  the  Markov 
chain  model,  if  we  could  reduce  the  dimension  from  n  to  m  «  n,  then  it 
would  significantly  speed  up  the  both  pairing  procedure  (  0(m{m  -  l)/2) 
instead  of  0(n(n  -  l)/2)  for  checking  pairs)  and  the  actual  sampling  from 
the  models.  In  fact,  this  dimension  reduction  even  may  reduce  the  amount 
of  unwanted  noise  in  the  images.  We  will  compare  the  performance  of  the 
reduced  models  with  that  of  the  full  models. 
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7.5  Integrating  Multiple  Resolutions  and  Sensors 

We  should  incorporate  all  the  pertinent  data  available  to  infer  the  loading 
condition  that  generated  the  observed  fracture  surfaces.  So  far ,  we  used  only 
the  surface  elevation  data  measured  by  the  FRASTAscope,  with  the  fixed 
magnification.  We  should  be  able  to  get  more  information  about  the  perti¬ 
nent  features  by  looking  at  different  resolutions.  We  will  also  examine  SEM 
images,  which  are  more  commonly  available  to  many  fractographers.  This 
availability  implies  that  a  large  empirical  rules  and  diagnostics  are  accumu¬ 
lated  among  the  experienced  fractographers.  Our  goal  here  is  to  interview 
these  fractographers  to  find  out  what  features  they  are  looking  for  to  infer 
the  loading  conditions  and  history.  The  FRASTAscope  can  also  generate  the 
contrast/reflectivity  images  of  fracture  surfaces.  These  images  may  be  used 
to  link  the  surface  elevation  data  and  the  SEM  images. 

7.6  Summary  of  our  recommendations 

The  above  projects  may  lead  to  a  better  classification  performance  and  better 
understanding  of  the  underlying  physics  of  the  fatigue  fractures,  which  is  the 
key  of  the  entire  project.  We  also  believe  that  this  project  will  generate 
significant  by-product  that  will  be  useful  for  other  area  of  DARPA’s  interest. 
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During  testing  of  a  high  pressure  pump  in  the  space  shuttle  fuel  tank,  a  crack  about  6  inches 
long  and  0.290  inch  deep  was  found  in  a  0.070-inch  blend  radius  of  the  turbine  housing  after  32 
starts.  A  failure  analysis  team  from  Pratt  &  Whitney,  Rocketdyne,  and  NASA  Lewis  looking  into 
the  cause  and  evolution  of  this  crack  wanted  to  know  when  in  the  test  history  the  crack  initiated,  if 
the  crack  was  unstable  in  the  late  stages  of  development,  and  whether  propagation  resulted  from 
engine  start/stop  cycles  or  power  level  changes.  To  answer  these  questions,  they  hoped  to  identify 
and  interpret  features  on  the  crack  surfaces. 

The  housing,  which  functions  as  a  pressure  vessel  for  high  temperature  (up  to  500°F) 
hydrogen  and  steam,  is  made  of  a  wrought  nickel-base  alloy  with  a  maximum  average  grain 
diameter  of  9  pm.  Chemical  analysis  showed  a  higher  concentration  of  hydrogen  near  the  surface 
than  in  the  bulk  of  the  material,  and  metallographic  sections  through  the  crack  showed  higher 
hydrogen  concentration  along  the  crack  and  at  the  crack  tip.  Engineering  estimates  of  the  strains 
indicated  levels  sufficiently  beyond  yield  that  linear  elastic  fracture  mechanics  principles  would  not 
apply.  Several  elastic-plastic  models  predicted  the  maximum  stress  intensity  at  the  tip  of  the  crack 
when  it  was  0.075  inch  long  to  be  43  to  77  ksiVim  A  of  77  ksi^im  should  produce  crack  . 
growth  in  hydrogen. 

Fractographic  features  were  not  distinctive  enough  to  characterize  the  stress  state,  but 
macroscopically  there  were  large  charges  in  the  surface  texture.  The  team  therefore  sought  other 
methods  to  analyze  the  fracture  surface  topography.  A  small  section  of  one  fracture  surface 
(Section  6B  from  Unit  9-1,  shown  in  Figure  B-l)  was  sent  to  SRI  for  analysis. 

SRI  used  the  FRASTAscope  to  produce  topographic  maps  of  29  areas  of  the  fracture 
surface  along  the  line  A-A.  Each  area  was  139.8  pm  in  the  direction  of  crack  growth  (i.e.,  in  the 
through-thickness  direction)  and  177.0  pm  in  the  perpendicular  direction.  This  area  was  divided 
into  a  600  by  800  pixel  array  and,  hence,  the  data  spacing  was  0.233  pm  in  both  directions. 

Figure  B-2  shows  contrast  images  and  gray-scale  topography  images  of  selected  areas. 
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Figure  B-1 .  Fracture  surface  from  a  Space  Shuttle  turbine  housing.  The  topographies  of  areas  along 
the  line  A-A  were  characterized  and  analyzed  by  FFT. 
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(a)  Contrast  image  of  the  area  at  69.3  pm  (b)  Gray-scale  topography  image  of  the 
from  the  inner  edge  area  at  69.3  pm  from  the  inner  edge 


(c)  Contrast  image  of  the  area  at  2379  pm  (d)  Gray-scale  topography  image  of  the 
from  the  inner  edge  area  at  2379  pm  from  the  inner  edge 
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Figure  B-2.  A  series  of  contrast  and  gray-scale  topography  images  of  selected  areas 
along  the  line  A-A  in  Figure  B-1 . 
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(e)  Contrast  image  of  the  area  at  2469  pm  (f)  Gray-scale  topography  image  of  the 
from  the  inner  edge  area  at  2469  pm  from  the  inner  edge 


(g)  Contrast  image  of  the  area  at  2869  pm  (h)  Gray-scale  topography  image  of  the 
from  the  inner  edge  area  at  2869  pm  from  the  inner  edge 
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Figure  B-2.  A  series  of  contrast  and  gray-scale  topography  images  of  selected  areas 
along  the  line  A-A  in  Figure  B-1 .  (Continued) 
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(i)  Contrast  image  of  the  area  at  4269  pm  (j)  Gray-scale  topography  image  of  the 
from  the  inner  edge  area  at  4269  pm  from  the  inner  edge 


(k)  Contrast  image  of  the  area  at  4469  pm  (I)  Gray-scale  topography  image  of  the 
from  the  inner  edge  area  at  4469  pm  from  the  inner  edge 
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Figure  B-2.  A  series  of  contrast  and  gray-scale  topography  images  of  selected  areas 
along  the  line  A-A  in  Figure  B-1 .  (Continued) 
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(m)  Contrast  image  of  the  area  at  4659  pm  (n)  Gray-scale  topography  image  of  the 
from  the  inner  edge  area  at  4659  pm  from  the  inner  edge 
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Figure  B-2.  A  series  of  contrast  and  gray-scale  topography  images  of  selected  areas 
along  the  line  A-A  in  Figure  B-1 .  (Concluded) 
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A  fast  Fourier  transform  analysis  was  applied  to  the  elevation  files,  and  the  averaged 
elevation  power  spectrum  density  (EPSD)  distribution  was  computed  for  each  area.  The  results  are 
plotted  in  Figures  B-3(a),  (b),  and  (c),  where  the  parameter  is  the  distance  from  the  centroid  of  the 
area  to  the  edge  of  the  housing.  Between  the  wavelengths  of  0.5  and  10  |im,  the  EPSD  curves  are 
roughly  parallel  to  one  another  and  slighdy  concave.  To  examine  the  variation  of  surface 
roughness  with  crack  depth,  we  plotted  the  value  of  the  ESPD  at  the  arbitrarily  chosen  wavelength 
of  5  pm  versus  distance  from  the  housing  edge,  Figure  B-4(a).  A  plot  of  the  square  root  of  the 
EPSD  versus  crack  depth  is  shown  in  Figure  B-4(b).  These  figures  show  two  large  peaks  in 
fracture  surface  roughness. 

In  light  of  the  results  presented  in  the  body  of  this  report,  these  peaks  suggest  loading 
excursions  at  crack  depths  of  2.5  mm  and  4.5  mm.  Because  reference  data  for  this  nickel-base 
alloy  have  not  been  generated,  we  are  unable  to  convert  the  EPSD  values  to  quantitative  estimates 
of  K,^  or  AK.  However,  the  slight  concavity  in  the  EPSD  curves  is  reminiscent  of  the  shape  of 
the  curves  for  the  titanium  alloy  reported  in  the  body  of  this  report.  Therefore,  one  might  speculate 
that  the  R  value  is  closer  to  0.1  than  to  0.8. 

Several  smaller  peaks  are  also  evident  in  Figure  B-4  and  suggest  smaller  overloads.  A 
similar  analysis  performed  on  center-cracked  panels  of  aluminum  alloy  sheet  showed  similar  EPSD 
peaks  at  crack  lengths  where  overloads  were  applied.9  Greater  peak  heights  were  obtained  with 
larger  overloads. 

It  would  be  of  great  interest  to  compare  the  results  of  Figure  B-4  with  the  run  time  history 
(starts/stops  and  power  changes  versus  time)  of  the  engine. 
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(c)  Elevation  Power  Spedum  Density  Plot  for  the  Last  Eight  Areas  along  the  Line  A-A 
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Figure  B-3.  Elevation  power  spectrum  density  plots  of  the  areas  along  the  line  A-A.  (concluded) 
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(a)  EPSD  values  as  a  function  of  crack  depth 


(b)  SQRT(EPSD)  values  as  a  function  of  crack  depth 
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Figure  B-4.  Variation  of  EPSD  and  SQRT(EPSD)  as  a  function  of  crack  depth. 
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